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CHAPTER  1  -  INTRODUCTION 

1.0  Problem  Formulation 

Adaptive  methods  of  beamforming  that  yield  high  resolution 
beam  outputs  require  a  much  greater  numerical  intensity  than 
conventional  Fourier  techniques  of  beamforming.  Reduction  of  the 
adaptive  dimension  of  these  data  dependent  beamformers  can 
alleviate  the  numerical  intensity  problem  with  a  negligible  loss  in 
performance. 

The  adaptive  dimension  reduction  may  be  accomplished  by  a 
matrix  premultiplication  of  the  array  sensor  data.  The  matrix 
premultiplication,  known  as  matrix  preprocessing,  transforms  the  N 
dimension  array  data  vector  into  an  M  dimensional  space  where  M  < 
N.  The  M  dimensional  data  is  then  beamformed  with  a  minimum 
variance  distortionless  response  (MVDR)  adaptive  beamformer. 

Another  method  of  reducing  the  adaptive  dimension  is  by  the 
approximation  of  the  covariance  or  cross  spectral  density  matrix 
(CSDM)  by  its  dominant  or  strongest  eigenvectors,  v  This  technique, 
developed  by  Owsley  in  [1]  and  [10],  is  known  as  the  enhanced 
minimum  variance  beamformer  (EMVDR). 

v 

The  intent  of  this  thesis  is  to  present  and  compare  the  matrix 
preprocessor  configuration  and  the  enhanced  minimum  variance 
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reduced  dimension  adaptive  beamformers.  The  subarray  and  beam 
space  matrix  preprocessors  will  be  covered.  The  comparison  will  be 
against  a  full  array  MVDR  beamformer  and  a  conventional 
beamformer. 

A  general  introduction  to  the  concept,  application  and 
implementation  of  beamforming  is  covered  in  Chapter  2. 
Conventional  and  adaptive  methods  of  beamforming  are  discussed 
and  the  reduced  dimension  adaptive  beamformer  is  proposed. 

The  models  and  assumptions  used  in  the  development  of  the 
beamformers  are  presented  in  Section  3.1.  Conventional 
beamforming,  the  beamformer  beampattern  and  beam  response  are 
covered  in  Section  3.2.  The  matrix  preprocessed  MVDR  algorithm 
development,  subarray  and  beam  space  preprocessors  and  a 
discussion  of  wavenumber  analysis  are  covered  in  Section  3.3.  The 
enhanced  minimum  variance  beamformer  is  developed  in  Section 
3.4. 


The  comparisons  will  first  be  in  the  form  of  a  metric  introduced 
by  Owsley  in  [4]  and  [10]  called  the  Array  Gain  Improvement  (AGI) 
in  Section  4.1,  then  beampatterns  in  Section  4.2,  followed  by  beam 
responses  in  Section  4.3  and  numerical  intensity  in  Section  4.4. 
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The  significant  contributions  of  this  thesis  include  the 
exploration  of  the  subarray  matrix  preprocessor  and  the  stochastic 
evaluation  of  the  beamformers. 

1.1  Current  Technology 

Adaptive  algorithms  have  origins  dating  to  the  1960’s.  The 
numerical  intensity  required  of  the  implementation  has  always  been 
a  hurdle  that  required  surmounting.  Alternatively  the  use  of  sub- 
optimal  techniques  has  been  proposed  in  [11].  With  the  onslaught  of 
VLSI  technology,  realization  of  broadband  full  array  adaptive 
beamformers  has  become  a  reality.  For  a  fixed  amount  of  processing 
power  a  trade-off  exists  between  the  size  of  the  array,  the  number  of 
beams  and  the  number  of  frequencies  the  beamformer  operates  on. 

The  MVDR  algorithm  was  first  developed  with  a  soft 
distortionless  response  constraint  by  Owsley  [17]  in  1969  and  with  a 
hard  constraint  as  shown  by  Frost  [18]  in  1972. 

The  reduced  dimension  adaptive  beamformers  to  be  discussed 
have  been  proposed  by  various  people.  The  subarray  preprocessor 
was  alluded  to  by  Owsley  in  [2],  the  beam  space  preprocessor  by 
Gray  in  [3]  and  the  enhanced  minimum  variance  beamformer  was 
developed  by  Owsley  in  [1]  and  [10].  Several  others  have  worked 
with  the  reduced  adaptive  dimension  idea  as  comprehensively 
shown  in  [12]  by  Van  Veen. 
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The  implementation  of  the  subject  beamformers  requires  two 
specialized  algorithms.  The  first  is  used  to  realize  the  covariance 
matrix  inversion  for  the  MVDR  beamformer  via  the  Cholesky  square 
root  factor  and  is  found  in  [2]  and  [5].  The  second  is  an  eigenvector 
and  eigenvalue  estimator  used  to  compute  the  eigen-decomposition 
of  the  covariance  matrix  for  the  EMVDR  beamformer  and  is 
developed  by  Yang  and  Kaveh  in  [7].  The  Cholesky  square  root  factor 
CSDM  exponential  update  and  inverse  and  the  eigenvector  algorithm 
by  Yang  and  Kaveh  are  respectively  shown  in  Appendices  A  and  B. 


CHAPTER  2  -  BEAMFORMING  BACKGROUND 


2.0  Introduction 

The  intent  of  this  section  is  to  introduce  the  reader  to  the 

concept  of  beamforming,  its  uses,  implementations  and  structures. 
The  properties  of  conventional  and  adaptive  beamforming  are 

discussed  and  the  concept  of  reduced  dimension  adaptive 
beamforming  is  proposed. 

2.1  Array  Signal  Processing 

Array  signal  processing  consists  of  the  processing  methods 
used  to  analyze  the  signals  received  from  an  array  or  group  of 
transducers.  The  signals  arriving  at  the  array  are  traveling  waves  of 
differing  frequencies  and  arrival  angles.  Applications  of  array  signal 
processing  may  be  found  in  seismology,  radar,  sonar,  and  other  fields 
as  seen  in  [1].  Examples  of  sensor  arrays  include  sensors  measuring 
seismic  activity  in  the  earth,  radar  antennae  and  underwater 
transducers  sensing  pressure  waves.  The  signals  received  by  the 

sensors  are  used  to  estimate  parameters  associated  with  the  sources 

causing  the  traveling  waves.  The  parameters  may  include 
frequencies,  amplitudes,  relative  phase  delays,  angles  of  arrival, 
source  location,  etc.  Beamforming  is  an  array  signal  processing 
technique  that  aids  in  the  estimation  of  these  parameters  and  in  the 
detection  of  the  traveling  waves. 
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2.2  Beamforming 

Beamforming  combines  the  data  from  the  array  sensors  into  a 
single  "beam"  output.  Several  beams  are  computed  so  that  they 
"look"  in  directions  that  span  the  possible  arrival  angles  of  an 
impinging  traveling  wave.  The  directional  aspect  of  beamforming  is 
accomplished  by  the  attenuation  and  phasing  or  delay  of  signals 
coming  from  different  directions.  This  process  may  be  related  to  a 
bandpass  finite  impulse  response  (FIR)  filter  operating  on  time  series 
data.  The  transfer  function  of  the  filter  is  similar  to  the 
beamformer's  beampattem  where  different  directions  are  attenuated 
instead  of  frequencies.  The  beampattern  is  discussed  more 
thoroughly  in  Section  3.2.  A  beamformer  is  thus  also  known  as  a 
spatial  filter.  The  relations  between  spectral  estimation  and 
beamforming  will  be  exploited  when  possible  to  assist  the  reader  in 
relating  the  text  to  a  familiar  topic. 

Other  characteristics  of  beamforming  include  the  increase  in 
signal  to  noise  ratio  (SNR)  at  the  output  of  a  beam  as  compared  to  the 
SNR  at  a  single  sensor.  The  ratio  of  the  beam  output  SNR  to  the 
sensor  level  SNR  is  known  as  the  array  gain  [8].  As  in  frequency 
estimation  the  angle  of  arrival  or  spatial  frequency  of  a  traveling 
wave  may  be  found  by  forming  several  beams  and  estimating  the 
spatial  spectrum. 
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2.3  Time  and  Frequency  Domain  Beamformers 

Beamforming  may  be  accomplished  in  either  the  time  domain 
where  the  sensor  outputs  are  delayed,  scaled  and  summed  or  in  the 
frequency  domain  where  the  delay  operation  translates,  via  a  Fourier 
transform,  into  a  phase  shift.  The  phase  shift  of  the  frequency 
domain  beamformer  is  accomplished  by  the  multiplication  of  a 
frequency  dependent  complex  exponential. 

The  time  domain  beamformer  operates  on  all  frequencies  and 
is  thus  wideband.  A  frequency  domain  beamformer  first  converts 
each  sensor  output  time  signal  into  the  frequency  domain  by  a  Fast 
Fourier  Transform  (FFT).  The  sensor  outputs  at  discrete  frequencies 
bins  are  then  scaled,  phase  shifted  and  summed  to  form  a  single 
frequency  beam  output.  The  beamforming  is  thus  performed  over  a 
frequency  band  as  wide  as  an  FFT  bin.  If  several  frequency  bins  are 
beamformed  and  converted  back  into  the  time  domain  by  an  inverse 
FFT  (IFFT),  then  time  domain  beam  outputs  are  available.  Figure  2.1 
describes  the  frequency  domain  beamformer  which  will  be  the 
implementation  discussed  in  this  paper. 

2.4  Sensor  Array  Geometry 

Without  a  loss  of  generality  we  will  restrict  the  configuration  of 
the  sensor  arrays  used  as  examples  to  be  a  set  of  equally  spaced 
sensors  placed  in  a  straight  line.  This  may  be  directly  related  to  the 
uniform  time  sampling  performed  on  a  voltage  or  current  waveform. 
In  this  case  the  waveform  being  sampled  is  the  signal  field  impinging 
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Figure  2.1  Frequency  domain  beamformer. 
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on  the  sensor  array.  The  sources  generating  the  traveling  waves  are 
assumed  to  be  point  sources  and  the  propagation  of  the  waves 
istaken  to  be  spherical.  Figure  2.2  depicts  a  line  array  of  sensors 
with  a  point  source  emitting  a  wave  impinging  on  the  array.  The 
perpendicular  direction  to  the  line  array  is  known  as  "broadside". 
From  broadside  traveling  waves  may  arrive  at  the  array  in  angles 
ranging  from  -90  to  +90  degrees. 

Source 


Figure  2.2  Line  array  with  point  source  at  angle  theta. 

2.5  Conventional  Beamforming 

Conventional  beamforming  (CBF)  is  beamforming  with  the  use 
of  fixed  amplitude  scalings  and  phase  shifts.  This  provides  a 
constant  beampattem  that  depends  on  the  scalings  used  and  the  total 
length  or  aperture  of  the  array.  This  is  the  most  widely  accepted 
type  of  beamforming.  The  spectral  estimation  analog  of  CBF  is  the 
discrete  fourier  transform  (DFT)  or  periodogram  method  with 
amplitude  shading  of  the  data.  Marple  calls  this  the  Welch 
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periodogram  in  [9].  Figure  2.3  shows  a  typical  beam  pattern  plot. 
The  pointing  direction  of  the  beam,  commonly  called  the  main 
response  axis  (MRA)  or  the  angle  the  beam  is  steered  to,  is  the  angle 
at  which  an  incoming  signal  receives  no  degradation. 


Figure  2.3  Conventional  beamformer  beampattem. 

2.6  Adaptive  Beamforming 

As  in  spectral  estimation  there  are  several  techniques  of 
beamforming  that  are  dependent  on  the  data.  These  algorithms  are 
known  as  adaptive  beamforming  (ABF).  Typically  the  ABF 
techniques  are  characterized  by  high  resolution  beams  and  often  an 
optimality  criterion.  Adaptive  beamformers  find  their  origin  in  the 
1960's  with  the  work  of  people  such  as  Applebaum,  Byrn,  Griffiths, 
et  al. 


The  adaptive  technique  used  in  this  paper  is  in  conjunction 
with  the  Minimum  Variance  Distortionless  Response  (MVDR) 
beamformer.  This  is  an  optimal  technique  that  minimizes  the  noise 


variance  at  the  output  of  a  beam  while  constraining  the  response  of 
the  beamformer  to  a  source  impinging  on  the  beam's  MRA  to  be 
distortionless,  i.e.  to  have  a  frequency  response  at  the  MRA  with 
constant  magnitude  and  zero  phase.  It  will  be  shown  that  the 
minimum  beam  output  noise  variance  condition  also  satisfies  the 
maximum  array  gain,  maximum  beam  output  SNR  and  minimum 
beam  output  power  conditions. 

Adaptive  beamforming  provides  an  increase  in  the  resolution 
of  a  beam  over  a  conventional  beamformer  under  certain  conditions. 
The  increased  resolution  assists  in  the  estimation  of  the  direction  of 
arrival  of  a  signal  impinging  on  the  array  of  sensors. 

In  order  to  implement  the  optimal  beamforming  algorithm 

some  information  about  the  traveling  waves  impinging  on  the  array 

is  necessary.  The  covariance  matrix  contains  that  information  and  if 
% 

it  is  not  known  a  priori  it  must  be  estimated.  The  data  dependent 
characteristic  of  adaptive  beamforming  is  in  the  estimation  of  the  N- 
by-N  (for  an  array  of  N  sensors)  covariance  matrix.  The  covariance 
matrix  is  covered  in  detail  in  Section  3.1.  In  reference  [2]  Owsley 
states  that  the  variance  of  the  exponentially  averaged  covariance 
matrix  estimator  is  inversely  proportional  to  S-N+l  where  S  is  the 
number  of  independent  updates  made  to  the  estimate  and  N  is  the 
number  of  sensors  (it  is  assumed  that  S  is  greater  than  N). 
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The  MVDR  ABF  algorithm  requires  the  inverse  of  the 
covariance  matrix  post  multiplied  by  a  known  vector.'  As  seen  in  [2] 
this  requires  on  the  order  of  N2  numerical  operations  utilizing  a 
Cholesky  square  root  factorization  of  the  covariance  matrix.  A 
numerical  operation  is  considered  to  be  a  complex  multiply  and  add. 
Comparatively  a  CBF  only  requires  N  numerical  operations  to 
compute  a  single  beam  output. 

The  costs  of  using  ABF  thus  include  a  noise  variance  on  each 
beam  output  due  to  the  uncertainty  of  the  estimated  covariance 
matrix  and  the  increased  number  of  numerical  operations  required. 

2.7  Reduced  Dimension  Adaptive  Beamforming 

There  have  been  several  studies  [1,  2,  3,  4  and  12]  concerning 
methods  of  reducing  the  dimension  of  the  adaptive  process  in  ABF, 
resulting  in  a  reduction  in  the  numerical  intensity  and  possibly  the 
estimation  noise  variance.  Three  of  the  more  promising  algorithms 
have  been  chosen  to  be  discussed  here:  The  subarray  preprocessor 
introduced  by  Owsley  in  [2],  the  beam  space  preprocessor  presented 
by  Gray  in  [3],  and  finally  the  enhanced  minimum  variance 
beamformer  developed  by  Owsley  in  [1]  and  [10]. 

The  subarray  and  beam  space  methods  are  both  characterized 
by  a  matrix  preprocessor  operating  on  the  array  sensor  data 
reducing  the  data  vector  size  and  the  adaptive  dimension  from  N  to 
M  as  shown  in  [4].  An  adaptive  MVDR  beamforming  process  is 
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performed  on  the  M  dimensional  output  of  the  matrix  preprocessor. 
A  reduction  in  the  number  of  numerical  operations  is  achieved 
because  of  the  reduced  dimension  of  the  adaptive  process. 

The  enhanced  minimum  variance  beamformer  (EMVDR)  uses 
an  estimate  of  the  covariance  matrix  based  on  its  dominant  modes  to 
compute  the  beamformer  weights  and  delays.  Only  the  D  dominant 
eigenvectors  must  be  estimated  as  compared  to  the  full  covariance 
matrix  where  D  is  expected  to  be  much  less  than  N  the  full 
dimension.  The  computational  savings  is  also  evident  in  the 
eigenvector  estimation  algorithm. 

The  purpose  of  this  thesis  is  to  develop  the  theory  behind  the 
above  three  reduced  dimension  adaptive  beamforming  techniques 
and  compare  them  to  conventional  beamforming  and  to  a  full  array 
adaptive  MVDR  beamformer. 


CHAPTER  3  -  THEORY 


3.0  Introduction 

This  section  is  intended  to  describe  the  models  and 
assumptions  used  in  the  development  of  the  beamformers  including 
the  introduction  of  the  "steering  vector"  and  the  covariance  matrix. 
Conventional  beamforming  is  discussed  along  with  a  beamformer's 
beampattem  and  beam  response.  The  MVDR  algorithm  with  a  matrix 
preprocessor  is  developed  followed  by  an  in  depth  discussion  of  the 
subarray  preprocessor.  The  beam  space  preprocessor  is  presented 
and  finally  the  enhanced  minimum  variance  beamformer  is  derived. 

The  variable  notation  used  throughout  the  text  has  the 
following  explanation:  scalars  are  generally  lower  case  non-bold 
characters,  vectors  are  lower  case  bold,  matrices  are  upper  case 
bold,  and  Fourier  transforms  may  be  either  upper  or  lower  case  of 
the  time  domain  signal  depending  on  the  situation. 

3.1  Models  and  Assumptions 

As  mentioned  in  Section  2.4  and  described  in  Figure  2.2  the 
example  used  to  develop  the  beamformers  will  be  an  array  of 
sensors  configured  in  a  straight  line  with  equal  spacing  between  each 
sensor.  The  sensor  is  a  transducer  that  converts  a  physical  force  into 
a  voltage  signal.  The  voltage  signal  is  then  sampled  in  time.  The 
resulting  discrete  voltage  signal  is  converted  into  the  frequency 
domain  by  the  use  of  an  FFT.  Figure  3.1  depicts  this  process.  A 
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vector  is  formed  by  taking  the  same  frequency  bin  (FFT  output)  from 
each  sensor  signal.  This  is  known  as  the  frequency  data  vector.  For 
wideband  processing  this  must  be  done  for  several  frequency  bins. 


Traveling  waves 
impinging  on 
array 


Time  Fast  Fourier 


Sensor  domain  Transform 


Figure  3.1  Pre-beamformer  signal  path. 


Figure  2.2  illustrates  a  point  source  producing  waves 
propagating  spherically.  It  is  assumed  that  all  point  sources  are  in 
the  far  field.  The  far  field  or  plane  wave  assumption  is  that  the 
length  of  the  array  of  sensors  is  small  compared  to  the  distance 
between  the  source  and  the  array.  If  this  is  so  the  wavefront,  as  it 
arrives  at  the  array,  is  very  near  to  being  planar.  Figure  3.2  depicts 
a  plane  wave  impinging  on  a  line  array  from  the  angle  theta.  The 
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resulting  time  delay  between  the  first  and  the  ith  sensor  is  seen  to 
be, 

dsin(e)(i  -  1)  „  . 


where  c  is  the  speed  of  of  the  traveling  wave  and  d  the  distance 
between  adjacent  sensors. 

dsin(  0) 

x  = - 

c 


Adjacent 
sensors 

Figure  3.2  Delay  between  adjacent  sensors. 

Grouping  all  of  the  sensor  outputs  due  to  the  signal  into  a 
vector  and  incorporating  the  delays,  Tj,  yields, 

a(t)  =  [a(t-xi)  a(t-X2)  -  a(t-xN)]T  , 


where  a(t)  is  the  amplitude  of  the  received  waveform  at  the  first  or 
reference  sensor  and  the  superscript  "T"  denotes  the  transpose 
operation.  Computing  the  Fourier  transform  of  the  time  domain  data 
vector  yields. 
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+®° 

a(f)  =  / a(t)  e-j2*ft  dt 

-oo 

a(f)  =  A(f)  [  e-i2n{x2  -  e-j2*riN]T 

a(f)  =  A(f)  d(f,0)  3-2 

where  A(f)  is  the  Fourier  transform  of  a(t).  It  is  assumed  that  A(f) 
has  zero  mean,  E[A(f)J=0,  and  variance  E[A(f)A*(f)]  =  S,  where  S  is 
frequency  dependent. 

The  vector  d(f,0)  is  dependent  on  the  frequency  and  the  angle 
of  arrival  of  the  plane  wave,  the  array  sensor  spacing  and  the  speed 

of  the  traveling  wave.  This  vector  is  known  as  the  steering  vector 

and  will  be  seen  to  play  an  important  role  in  beamforming.  The 
steering  vector  is  developed  here  for  the  simple  line  array  case. 

More  complex  array  configurations  only  require  that  the  relative 
delay  between  a  reference  sensor  and  every  other  sensor  be 

calculable  or  measurable. 

In  reality  the  sensors  receive  noise  as  well  as  signals  from  the 
ambient  background  and  electronic  noise  from  amplifiers  and  other 
hardware.  This  noise  is  modelled  as  being  uncorrelated  from  sensor 
to  sensor  and  additive  to  the  received  voltage  signal,  i.e. 


x(f)  =  a(f)  +  n(f) 


3-3 


1  8 

where  n(f)  is  the  Fourier  transform  of  the  received  noise.  The 
frequency  noise  vector  is  characterized  by  having  a  complex  zero 

mean  (E[n(f)]=n(f)=0)  and  a  covariance  matrix  defined  by: 

Q  =  E[(n  -  n)(n  -  n)H]  =  E[nnH]  =  a2  In 

where  the  superscript  "H"  denotes  a  complex  conjugate  transpose 
operation  and  In  is  an  N-by-N  identity  matrix.  For  ease  of  notation 
the  dependence  on  frequency  will  be  suppressed  for  some  variables. 

The  covariance  matrix  mentioned  in  Chapters  1  and  2  is  the 
expected  value  of  the  outer  product  of  the  frequency  data  vector 
with  itself, 

R  =  E[x  xH] 

=  E[(a  +  n)(a  +  n)«] 

=  E[aaH  +  anH  +  naH  +  nnHl 

Since  a  and  n  are  assumed  statistically  independent,  the  expected 
value  of  their  outer  product  is  zero  yielding, 

=  E[aaH  +  nnH) 

=  E[A(f)A*(f)  ddH]  +  a2  IN 
=  E[A(f)A*(f)]  ddH  +  a2  IN 
=  S  ddH  +  a2  jN  . 


3-4 
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Now  let  there  be  an  arbitrary  number  of  independent  or 

uncorrelated  plane  wave  sources  impinging  on  the  array  from 

directions  61,  02,  0p  with  the  associated  steering  vectors  d(0i), 
d(02),  •••  d(0p)  and  source  powers  (variances)  Si,  S2,  —  Sp.  Due  to 

linearity  and  superposition,  the  terms  resulting  from  the  plane  wave 

sources  add  to  form  a  covariance  matrix  of  the  following  form: 

R  =  E[x  xH] 

=  Si  didiH  +  S2  d2d2H  +  —  +  Sp  dpdpH  +  In 

=  [di  d2  —  dp]  S  [d  1  d2  —  dp]H  +  In 

=  D  S  DH  +  a2  IN  3-5 

where  D  =  [di  d2  •  dp],  and  S  =  diag[Si  S2  •••  Sp]. 

The  covariance  matrix,  described  by  the  expected  value  of  the 
outer  product  of  the  frequency  data  vector  with  itself,  is  also  known 
as  the  cross  spectral  density  matrix  (CSDM)  for  a  particular 
frequency. 

3.2  Conventional  Beamforming 

Conventional  beamforming  is  performed  with  fixed  amplitude 
scalings  and  phase  shifts  that  correspond  to  fixed  look  directions. 
The  spectrum  analysis  analog  is  the  bandpass  FIR  filter.  The  weights 
applied  to  the  sensors  are  similar  to  the  impulse  response  or 
coefficients  of  the  FIR  filter.  The  beam  output  is  the  sum  of  the 
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products  of  the  complex  sensors  weights  and  the  sensor  data,  i.e.  the 
discrete  convolution  of  the  two. 

The  spectral  analysis  technique  known  as  the  Welch 
periodogram  [9]  is  performed  using  the  DFT  which  is  simply  a  bank 
of  FIR  bandpass  filters.  Thus  it  is  seen  that  the  FIR  spatial  filter  (the 
CBF)  applied  to  the  sensor  array  data  is  simply  estimating  the  spatial 
frequency  spectrum  (at  a  single  spatial  frequency)  via  a  Welch 
periodogram  technique.  It  is  convenient  to  think  of  a  CBF  as  a  spatial 
bandpass  FIR  filter.  In  general  several  beams  are  computed  to  form 
an  estimate  of  the  full  spatial  spectrum.  The  spectrum  analysis 
analog  is  a  bank  of  FIR  filters  estimating  several  frequencies 
(probably  implemented  as  an  FFT). 

In  spectral  estimation  the  resolution  achieved  by  digital 
Fourier  techniques  such  as  periodograms  is  inversely  proportional  to 
the  amount  of  time  data  that  is  processed.  A  similar  limit  is  imposed 
in  conventional  beamforming  where  the  amount  of  data  processed  is 
defined  by  the  aperture  or  extent  of  the  array.  The  adaptive 
algorithms,  as  in  spectral  estimation,  are  not  restricted  by  this  limit. 

The  beampattern  and  beam  response  of  a  beamformer  are 
commonly  misunderstood.  The  beampattern  is  similar  to  the 
transfer  function  of  an  FIR  filter.  It  is  defined  as  being  the  expected 
response  of  a  beamformer  steered  to  a  specific  direction  (i.e.  a  fixed 
set  of  beamforming  weights)  due  to  a  single  unit  variance,  zero  mean. 


plane  wave  source  being  swept  across  all  potential  angles  of  arrival 
[5].  The  transfer  function  of  an  FIR  filter  may  be  measured  similarly 
by  sweeping  a  sine  wave  across  all  possible  frequencies.  The 
information  the  beampattem  yields  is  the  attenuation  that  a  plane 
wave  signal  may  expect  to  receive  dependent  on  its  angle  of  arrival 
when  the  beamformer  is  steered  to  a  fixed  angle.  As  previously 
mentioned  the  beampattern  provides  the  "impulse"  response  or 
transfer  function  of  the  beamformer  that  yields  the  output  or  beam 
response. 

The  frequency  domain  beamformer  is  represented  by  the  inner 
product  of  a  beamforming  vector  of  complex  weights,  w(0s),  with  the 
received  frequency  data  vector, 

y(6s)  =  w(0s)Hx(f)  . 

The  output  of  the  beamformer  steered  to  the  angle  0S  is  y(0s).  Note 
that  y  is  a  complex  scalar  value.  The  magnitude  squared,  or 
instantaneous  power,  of  the  beam  output  is  often  used  as  a 
comparison  metric  or  for  display  purposes. 

Mathematically  the  beampattern  is  developed  from  a  fixed  set 
of  weights,  w(0s),  where  0S  is  the  angle  the  beamformer  is  steered  to. 
The  weights  may  come  from  any  type  of  beamformer.  The  expected 
value  of  the  magnitude  squared  of  the  beamformer  output  is  the 
beampattern. 
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b(es,0)  =  E[iw(es)HX(f,e)i2]  3  -  6 

where  x(f,0)  represents  the  frequency  domain  data  vector  due  to  a 
unit  variance,  zero  mean,  plane  wave  emanating  from  the  angle  0. 
As  seen  in  Section  3.1  the  frequency  domain  data  vector  due  to  a 
plane  wave  source  is  modeled  as  a  random  complex  scalar  multiplied 
by  a  steering  vector  for  the  array  steered  to  the  direction  of  the 
source, 

x(f,0)  =  A(f)d(0) , 

where  the  complex  multiplier,  A(f),  has  zero  mean,  E[A(f)]  =  0,  and 
unit  variance  E[IA(f)l2]  =  1.  Now  simplifying  the  beampattem, 

b(0s,0)  =  E[w(0s)Hx(f,0)x(f,0)Hw(0s)] 

=  E[A(f)A*(f)w(0s)Hd(0)d(0)Hw(0s)] 

=  E[A(f)A*(f)]  w(0s)Hd(0)d(0)Hw(0s) 

=  E[IA(f)|2]  lw(0s)Hd(0)|2 

=  lw(0s)Hd(0)|2  3-7 

The  direction  of  the  unit  variance,  zero  mean,  plane  wave 
source,  0,  is  swept  through  all  possible  arrival  angles  (-90  to  +90 
degrees).  The  resulting  plot  of  b(0s,0)  against  the  source  arrival 
angle,  0,  is  the  beampattem  of  the  beamformer  steered  to  the  angle 
0S. 
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A  conventional  beamformer  delays,  scales  and  sums  the  sensor 
outputs.  The  delays  are  determined  by  the  direction  the 
beamformer  is  steered  to  and  as  seen  in  Section  3.1  are  a  function  of 
the  array  spacing,  speed  of  the  traveling  wave  and  steering  angle 
(for  the  line  array  example).  The  shadings  utilized  are  similar  to 
those  used  in  digital  filter  design.  The  uniform  and  Chebyshev  (with 

4 

30dB  sidelobes)  shadings  will  be  utilized  here.  Figures  3.3  and  3.4 
depict  the  beampattern  of  a  uniformly  and  Chebyshev  shaded 
conventional  beamformer  steered  to  an  angle  of  30  degrees. 


-80  -60  -40  -20  0  20  40  60  80 

Angle  of  arrival  (degrees) 


Figure  3.3  Uniform  CBF  beampattern  steered  to  30  degrees. 

The  characteristics  of  the  Chebyshev  shading,  the  constant 
sidelobe  level,  lower  than  that  of  uniform  sensor  shading  at  the 
expense  of  mainlobe  width,  are  seen  in  Figure  3.4. 
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The  beamforming  vector  for  a  conventional  beamformer, 
steered  to  the  angle  6S,  as  represented  by  the  delays  and  shadings  is 

w(0s)  =  W  d(6s) 

where  W  is  a  diagonal  matrix  containing  the  shadings  and  d(6s)  is  an 
array  steering  vector  steered  to  the  angle  0S.  The  beampattem  for  a 
conventional  beamformer  is  thus  seen  to  be 

bcBF(0s,e>  =  ld(0s)HWd(0)l2  3  -  8 

where  W  is  real  and  diagonal. 


■80  -60  -40  -20  0  20  40  60  80 

Angle  of  arrival  (degrees) 


Figure  3.4  Chebyshev  CBF  beampattem  with  30  dB  sidelobes 

steered  to  30  degrees. 
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The  beam  response  of  a  beamformer  is  the  output  power  of  a 
fan  of  beams  spanning  the  potential  arrival  angles.  The  spectrum 
analysis  analog  is  the  estimate  of  the  power  in  several  frequency 
bins.  The  beam  response  is  computed  by  forming  several  beams  and 
taking  the  magnitude  squared  of  the  output.  The  average  or 
expected  response  includes  the  expectation  operator  as  follows, 

r(0)  =  E[Iw(0)Hx|2] 

=  E[w(0)HxxHw(0)] 

=  w(0)H  E[xxH]  w(0) 

=  w(0)HRW(0)  3-9 

where  x  is  the  received  frequency  data  vector,  R  is  the  CSDM  defined 
in  Section  3.1,  w(0)  is  the  beamforming  vector,  and  0  is  the  angle  to 
which  each  beam  is  steered  to  and  varies  from  -90  to  +90  degrees. 

The  average  beam  response  provides  the  response  of  the 
beamformer  to  a  particular  scenario  for  the  perfectly  known  CSDM. 
In  practice  the  CSDM  is  never  known  exactly  and  scenarios  are  not 
stationary  for  long  periods  of  time.  The  CSDM  is  therefore  estimated 
using  an  exponential  averager  (see  [2],  [5]  and  Appendix  A).  The 
beam  response  for  each  iteration  is  calculated  as  the  power  output  of 
the  beamformer, 


r(0)  =  !w(0)Hxl2. 
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Figures  3.5  and  3.6  show  average  beam  responses  for  a 
conventional  beamformer  with  uniform  and  Chebyshev  (with  30  dB 
sidelobes)  shading.  The  scenario  is  a  plane  wave  source  at  zero 
degrees  with  OdB  variance  and  at  30  degrees  with  -3dB  variance 
amidst  spatially  uncorrelated  background  noise  with  a  variance  of 
OdB  at  a  sensor.  There  are  eleven  sensors  in  the  array. 
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Angle  of  arrival  (degrees) 


Figure  3.5  CBF  beam  response  with  uniform  shading. 

The  average  beam  response  for  a  conventional  beamformer  is 
seen  to  be, 

tcbfO)  =  d(0)HW  R  W  d(6),  3-11 

by  substituting  in  for  the  conventional  beamforming  vector. 
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Figure  3.6  CBF  beam  response  with  30  dB  Chebyshev  shading. 

3.3  Matrix  Preprocessed  MVDR 

The  subarray  and  beam  space  preprocessors  may  be  described 
in  the  form  of  a  matrix  premultiplier  operating  on  the  received 
frequency  data  vector.  The  implementation  of  this  configuration  is 
depicted  in  Figure  3.7  where  the  beamformer  is  separated  into  a 
preprocessor  and  an  ABF  section.  Note  that  the  full  array  MVDR 
algorithm  is  a  special  case  where  the  matrix  preprocessor  is  an  N-by- 
N  identity  matrix. 


N  =  11 

Chebyshev  shading 
with  30dB  sidelobes 
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The  objective  of  the  matrix  preprocessor  is  to  reduce  the 
dimension  of  the  adaptive  part  of  the  beamformer  from  N  (the 
number  of  sensors)  to  M.  The  effect  of  the  reduction  lies  in  the  fact 
that  the  CSDM  that  needs  to  be  estimated  and  inverted  is  now  only 
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M-by-M  thus  reducing  the  total  number  of  numerical  operations 
required  to  implement  the  beamformer.  Sections  3.3.1  and  3.3.2 
cover  two  particular  methods  of  choosing  the  matrix  used  in  the 
preprocessor.  As  shown  in  [4]  by  Owsley  and  Abraham  several  types 
of  reduced  dimension  adaptive  beamformers  may  be  implemented  in 
the  matrix  preprocessor  structure. 


•Preprocessor 


H, 


y(k,f)  =  A  x(k,f) 


Frequency 
data  vector 


Adaptive 
Linear  Combiner 


r(k,f)  =  w(k,f,0)  y?k,f) 


r(k,f,0) 

Beam 

output 


Figure  3.7  Matrix  preprocessor  -  ABF  configuration. 


The  minimum  variance  distortionless  response  beamformer 
may  be  equivalently  derived  for  the  constrained  minimization  of  the 
variance  of  a  beam  output,  the  maximization  of  a  metric  known  as 
the  array  gain,  the  maximization  of  the  beam  output  SNR  or  the 
minimization  of  the  total  beam  output  power.  Beginning  with  the 
maximization  of  the  array  gain,  the  equivalency  will  be  shown.  First 
the  frequency  data  vector  is  acted  on  by  the  preprocessor  matrix 
yielding  a  reduced  dimension  data  vector. 
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y(6ps)  =  AH(0ps)  x  3-12 

(M  x  1)  (M  x  N)  (N  x  1) 

where  8ps  is  the  presteering  angle  and  A(8ps)  is  the  preprocessing 
matrix.  Several  presteering  angles  are  required  to  evaluate  all  the 
possible  arrival  angles.  The  presteering  angles  are  specific  to  the 
type  of  preprocessor  matrix  used  and  will  be  discussed  further  in 
Sections  3.3.1  and  3.3.2. 

The  adaptive  beamforming  will  be  performed  on  the  reduced 
dimension  data  vector  y(9ps)*  thus  the  beam  output  power  or 
response  of  a  beam  steered  to  angle  8  is, 

r(8)  =  E[lw(0)Hy(eps)|2]. 

Note  that  the  beam  is  pointed  to  the  angle  0  and  the  preprocessor 
matrix  is  set  to  the  presteer  angle  9ps.  For  convenience  most  of  the 
indices  denoting  these  angles  will  be  omitted. 

r(8)  =  E[wHyyHw] 

=  wH  E[yyH]  w 
=  wH  E[AHx  xHA]  w 
=  wHAH  E[x  xHl  A  w 
=  wHAH  R  A  w 
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In  Section  3.1  the  CSDM  was  modeled  as  the  sum  of  several  plane 
waves  and  spatially  uncorrelated  noise, 

R  =  Si  didiH  +  S2  d2d2H  +  —  +  Sp  dpdp1*  +  a2  In  . 

Let  the  plane  wave  denoted  by  angle  81  be  from  the  angle  that  the 
beamformer  is  steered  to,  0  (i.e.  0  =  0i).  Now  an  important  distinction 
between  signal  and  interference  is  made.  Any  signals  impinging  on 
the  array  from  a  direction  other  than  the  look  direction  or  MRA  of 
the  beamformer  are  considered  to  be  interferences.  The 
uncorrelated  background  noise  is  also  considered  to  be  interference. 
Thus  the  CSDM  may  be  broken  up  into  signal  and  noise  partitions  as 
follows, 

P  =  djdiH  =ddH 

Sn  Q  =  S2  d2d2H  +  —  +  Sp  dpdpH  +  a2  In 

where  P  is  the  signal  partition,  Q  is  the  noise  partition  and  Sn  is  the 
equivalent  variance  at  a  single  sensor  due  to  all  of  the  interfering 
signals.  Now  the  CSDM  may  be  written  as, 

R  =  SjP  +  SNQ.  3-13 

Substituting  this  back  into  the  beam  response  power  equation  yields 


r(0)  =  wHAH  (SiP  +  SnQ)  Aw 


3  1 


=  SiwHAHPAw  +  SnwhAhQA  w  3-14 

It  is  desired  to  maximize  the  array  gain  (AG).  The  array  gain  is 
defined  by  Burdic  in  [8]  as  the  ratio  of  the  SNR  at  the  output  of  the 
beamformer  to  the  SNR  at  a  single  sensor. 


_  SNRqut 
~  SNRin 


3-15 


The  SNR  at  a  single  sensor  is  simply 

sNRm=|i . 
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The  SNR  of  the  beam  response  output  is  seen  to  be  the  ratio  of  the 
signal  and  noise  partitions  of  the  above  equation  for  the  beam 
response  output  power. 


SiwHAHP  Aw 
SNR0ut-SnWhahqAw- 


3-17 


The  array  gain  is  then 


_  wHAHP A w 
"whAHQAw  ’ 


3-18 


The  constraint  applied  to  the  beamformer  keeps  a  source 
coming  from  the  angle  that  the  beam  is  pointed  to  distortionless. 
This  means  that  the  signal  must  receive  no  phase  change  and  only  a 
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specified  amount  of  gain.  Here  the  gain  is  chosen  as  N,  the  number  of 
sensors.  This  is  to  insure  that  the  ABF  array  gain  matches  the  CBF 
array  gain.  The  constraint  is  written  as. 


wHAHd  =  N.  3-19 

The  steering  vector  d  and  the  beamforming  vector  w  are  pointed  to 
the  direction  of  the  beam,  and  in  this  case  that  of  the  first  plane 
wave  source.  A  cost  function,  Ji,  is  generated  that  includes  a 
Lagrange  multiplier,  X,,  operating  on  the  constraint, 

wHAHP  A  or 

Jl  =wHAHQAw  +  2MwHAHd  -  N).  3-20 

The  beamforming  vector,  w,  that  maximizes  the  cost  function  also 
maximizes  the  array  gain  while  preserving  the  distortionless 

response  constraint.  In  general  the  minimization  or  maximization  of 
the  cost  function  ;equires  that  the  partial  derivative  be  taken  with 
respect  to  the  beamforming  vector,  w,  and  its  complex  conjugate,  w*, 
and  then  setting  both  equal  to  the  zero  vector.  When  the  cost 
function  is  real  it  is  sufficient  to  set  either  partial  derivative  to  the 
zero  vector.  The  cost  function  is  real  here  because  the  function  being 
extremized  is  real  and  the  constraint  is  forced  to  be  real.  The 

multiplication  of  the  constraint  by  2  doesn't  alter  the  extremization 
of  the  cost  function  for  if  the  constraint  is  satisfied  (wHAHd  -  N)  will 

equal  zero  ana  the  factor  2  provides  a  simpler  algebraic  equation 

than  would  otherwise  have  occurred. 
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The  maximization  of  the  beam  output  SNR  condition  is  easily 
seen  to  be  the  previous  equation  for  the  array  gain  multiplied  by  the 
constant  sensor  level  SNR  yielding  a  cost  function. 


h  =  SNRin 


wHAHPAw 

wHAHQAw 


+  2X(wHAHd  N). 
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The  minimization  of  the  beam  output  noise  variance  may  be 
seen  by  substituting  the  constraint  into  the  expanded  numerator  of 
the  equation  for  the  array  gain. 


J3  = 


wHAHddHA  w 
wHAHQAw 


+  2Jt(wHAHd  -  N) 


wHpQAw  +  2X<wHAHd  -  • 
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The  minimization  of  the  beam  output  noise  power  or  variance  is 
equivalent  to  the  maximization  of  its  inverse.  The  denominator  of 
the  above  cost  function  is  the  power  output  of  the  beam  due  to  the 
noise  or  interference  portion  of  the  signals  impinging  on  the  array. 


The  minimization  of  the  total  beam  output  power  may  be  seen 
by  rewriting  the  cost  function  J3  as  a  minimization  rather  than  an 
maximization, 


J4  =  wHAHQAw  +  2X(wHAHd  -  N). 
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Now  substituting  in  for  the  noise  partition,  Q,  and  the  constraint  on 
the  beamforming  vector,  w, 


J4  =  wHAHRAw  -  wHAHddHAw  +  2X(wHAHd  -  N) 

=  wHAHRAw  -  N2  +  2X(wHAHd  -  N).  3-23 


The  constants  in  J2,  J3  and  J4  due  to  the  application  of  the 
constraint  or  the  sensor  level  SNR  do  not  alter  the  resulting  optimal 
beamforming  vector,  w.  Note  that  all  cost  functions  constrain  the 
response  of  a  signal  from  the  direction  the  beam  is  pointed  in  to  be 
distortionless  with  a  gain  of  N.  The  functions  to  be  extremized  are 
listed  in  Table  3.1  for  comparison. 


Table  3.1  Functions  to  be  extremized  as  related 
to  various  optimal  conditions. 


Description 

Function 

Extremum 

Array  gain 

wHAHP  A  w 
wHAHQAw 

Maximize 

Beam  output 
SNR 

wHAHPAw 

SNRinwhahqaw 

Maximize 

Beam  output 
Noise  power 

wHAHQAw 

Minimize 

Beam  output 

wHAHRAw  -  N2 

Minimize 

power 

The  optimal  beamforming  vector  is  most  easily  found  by 
solving  the  constrained  minimization  of  the  total  beam  output  power 
in  the  cost  function  J4.  The  partial  derivative  of  the  cost  function 
with  respect  to  the  beamforming  vector  is. 


a 

—  J4  =  2AhRAw  +  2XAHd 
dw 


Now  setting  this  equal  to  zero  and  solving  for  the  beamforming 
vector,  w. 


w  =  -X(AHRA)**AHd  3-24 

This  equation  for  w  is  substituted  into  the  constraint  equation  which 
is  then  solved  for  the  Lagrange  multiplier  X, 


wHAHd  = N 

-X  dHA(AHRA)*iAHd  =  N 

x _ _ 

A-‘dHA(AHRA)-lAHd  ' 

Note  that  (AHRA)*1  is  complex  conjugate  symmetric.  Substituting 
this  back  into  the  equation  for  w  thus  yields, 

N(A»RA)-iAHd 

W  "dHA(AHRA)-iAHd 

N  RP-1AHd 

-dHARPlAHd  3-25 

RP  =  AHRA 
=  AHE[xxH]A 
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=  E[yyH] 


where  y  =  AHx  is  the  reduced  dimension  frequency  data  vector  and 
Rp  is  the  reduced  dimension  CSDM  defined  as  the  expected  value  of 
the  outer  product  of  the  reduced  dimension  data  vector  with  itself. 

The  reduced  dimension  CSDM,  Rp,  may  be  estimated  from  the 
reduced  dimension  frequency  data  vector,  y.  There  exist  several 
methods  of  estimating  the  CSDM  of  a  frequency  data  vector  [5].  The 
method  used  in  Section  4.3  for  the  stochastic  analysis  of  the 
beamformers  is  an  exponential  average  of  a  rank  one  update  of  the 
CSDM  [2]  &  [5J.  In  particular  the  Cholesky  square  root  factor  of  the 
CSDM  is  what  is  updated.  This  method  is  covered  in  detail  in  [5]  and 
in  Appendix  A. 

The  development  of  the  MVDR  beamformer  is  covered  in  detail 
in  [1]  and  [5].  The  matrix  preprocessor  development  is  a  simple 
extension  that  came  from  work  done  for  [4]. 

3.3.1  Subarray  Preprocessing 

Subarray  preprocessing  consists  of  forming  several  smaller 
overlapping  arrays  out  of  all  of  the  sensors  of  the  full  array.  The 
smaller  arrays,  known  as  subarrays,  are  then  beamformed  using 
conventional  techniques  to  form  several  beam  outputs.  A  beam 
output  from  each  subarray  (steered  to  the  same  direction)  is 
adaptively  combined  according  to  the  previous  section's  development 
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of  the  matrix  preprocessed  MVDR  beamformer.  The  concept  of 
forming  subarrays  prior  to  adaptive  beamforming  is  alluded  to  by 
Owsley  in  [2].  Figure  3.8  shows  how  an  array  may  be  partitioned 
into  subarrays,  each  with  a  beam  pointed  to  a  certain  direction. 


Subarrays 


Subarray  beams 


Figure  3.8  Full  array  partitioned  into  subarrays. 


The  preprocessing  matrix  is  comprised  of  conventional 
beamforming  vectors  placed  such  that  they  pick  off  the  sensors  in 
each  subarray  and  combine  them  to  form  the  beam  outputs.  A 
typical  subarray  preprocessor  matrix  will  have  the  form. 


A(0P)  = 


w(9p) 

0 


0 

w(0p) 

0 


0 

w(0p)_ 


3-26 
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where  w(0p)  is  a  conventional  beamforming  vector  that  is  the  length 
of  a  subarray,  0p  is  the  direction  that  the  subarray  beams  are  steered 
to,  and  0  is  a  variable  dimension  vector  of  zeros. 

A  particular  concern  in  this  method  of  preprocessing  lies  in  the 
choice  of  the  subarray  size,  shading  and  spacing  of  the-  phase  centers. 
The  energy  field  the  sensors  are  immersed  in  is  sampled  in  discrete 
locations.  This  sampling  process  is  similar  to  the  sampling  of  a  time 
waveform  where  aliasing  is  avoided  by  sampling  above  the  Nyquist 
sample  rate.  In  time  waveforms  the  signal  is  filtered  prior  to  the 
sampling  to  insure  that  any  energy  at  frequencies  that  would  alias 
during  the  sampling  is  attenuated  as  much  as  possible.  Spatial 
sampling  does  not  have  this  luxury  and  hence  aliasing  may  not  be 
avoided. 

Throughout  the  paper  the  spatial  analog  of  temporal  frequency, 
spatial  frequency,  has  been  alluded  to  as  an  indication  of  the  angle  of 
arrival  of  a  traveling  wave  impinging  on  an  array.  In  order  to 
adequately  discuss  the  aliasing  concerns  of  the  subarray  parameters 
more  depth  in  the  spatial  frequency  domain  is  required.  Section 
3. 3. 1.1  covers  the  required  information. 

3.3. 1.1  Wavenumber  Analysis 

The  spacing  of  the  sensors  in  an  array  is  chosen  such  that  the 
highest  "spatial  frequency",  known  as  wavenumber,  that  is  desired  to 
be  evaluated  is  not  aliased  at  the  highest  temporal  frequency  of 


concern.  The  spatial  analog  to  time  and  period  is  distance  and 
wavelength, 

X  =  j  ,  3-27 

where  c  is  the  speed  of  the  propagating  wave  and  f  is  the  frequency 
of  concern.  Spatial  frequency  is  measured  in  wavenumber  defined 
by, 

2rcf  2tc 

k  =  —  sin(0)  =  —  sin(0)  ,  3-28 

for  wave  propagation  along  a  single  axis  (i.e.  a  line  array).  This  is 
covered  by  Owsley  in  [5]  for  the  more  general  three  dimensional 
wave  propagation  as  well  as  the  single  axis  case. 

The  sampling  in  space,  as  in  time,  causes  the  spectrum  of  the 
sampled  signal  to  be  the  original  unsampled  spectrum  repeated  at 
every  multiple  of  the  spatial  sampling  rate,  ks.  The  sensors  in  the 
array  only  cover  a  finite  portion  of  the  waves  traveling  in  the 
medium.  This  is  modeled  as  an  infinite  function  windowed  to  a  finite 
aperture.  The  resulting  spectrum  is  thus  the  convolution  of  the 
transform  of  the  windowing  function  with  the  true  spectrum.  The 
sampling,  windowing  and  convolution  principles  can  be  found  in 
Oppenheim  and  Schafer  [14]. 
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Beamforming  the  frequency  data  vector  may  be  viewed  as  a 
method  of  forming  an  estimate  of  the  spatial  spectrum  at  a  single 
frequency  and  a  particular  wavenumber.  The  frequency  and  the 
angle  the  beam  is  steered  to  define  the  wavenumber.  The  steering 
vector  used  in  conventional  beamforming  is, 

d(f,9)  =  [  c>2*h2  ...  e-j2'f,N]T 

=  [  1  e-i2*Msi"«»/c  ...  e-j2jtfdsin(6)(N-l)/CiT 
=  [e-j2”Msin(8)i/cj  .  =  a  N., 

=  [e‘jbi] 

Looking  at  the  exponent  only, 

bj  =  27tfdsin(0)i/c 
2*f 

=  —  sin(0)  d  i 

=  kd i  3-29 

The  sensor  spacing  of  the  array  is  chosen  such  that  it  is  no 
greater  than  one-half  of  the  wavelength  of  the  highest  wavenumber 
of  concern.  This  is  related  to  the  frequency  that  the  array  is 
designed  for  and  the  speed  of  the  traveling  waves  in  the  medium, 

2  =  2fs  =  ks  • 


41 


Thus  substituting  this  back  into  the  above  equation  for  the  exponent 
yields. 


The  wavenumbers  to  be  evaluated  that  will  produce  unique  spatial 
spectrum  estimates  are  seen  to  lie  on  [-ks,  ks]  where. 


This  may  be  shown  by  evaluating  a  wavenumber  with  a  greater 
magnitude  than  ks, 

k'  =  k  +  ks, 

where  k  is  less  than  ks.  The  exponent  bi  is  then. 


bi  =^7™  ‘ 


Substituting  this  back  into  the  steering  vector  causes  the  "27ti"  term 
to  go  away  because  it  is  equivalent  to  a  multiplicative  one.  The  value 
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of  wavenumber  that  is  thus  evaulated  by  using  k',  greater  than  ks,  is 
-(ks-k)  which  is  within  the  range  [-ks,  ks]. 

The  less  than  half  wavelength  spatial  sampling  may  be  related 

T 

to  spectral  analysis  where  the  time  sampling  occurs  every  ^  seconds 

! 

or  less.  The  unique  frequencies  are  evaluated  over  [-fs,  fs],  where  fs 
2 

£  y  .  In  this  case  the  Nyquist  sampling  frequency,  or  the  frequency 
interval  at  which  the  repeated  spectra  occur  is  fN  =  2fs. 

Figure  3.9  shows  the  wavenumber  -  frequency  spectrum  and 
the  repeated  spectra  at  multiples  of  the  spatial  Nyquist  sampling 
rate,  2ks.  Note  that  for  a  fixed  travelling  wave  speed,  c,  all  the 
spatial  energy  falls  within  the  shaded  region  where  the  boundaries 
indicate  a  wave  impinging  on  the  array  at  +90  or  -90  der.-ees  from 
broadside. 

The  beamforming  process  estimates  the  wavenumber  spectrum 
inside  of  the  shaded  region  seen  in  Figure  3.9  at  a  single  frequency. 
The  angles  the  beams  are  steered  to  correspond  to  particular 
wavenumbers.  Wavenumber  is  related  linearly  to  frequency  as  seen 
in  the  equation  relating  the  two  for  a  fixed  angle  of  arrival. 

3. 3. 1.2  Subarrays 

The  conventional  beamforming  of  the  subarrays  may  be 
thought  of  as  applying  a  wavenumber  bandpass  finite  impulse 
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response  (FIR)  filter  to  the  wavenumber  spectrum.  The  output 
spectrum  is  bandlimited  in  wavenumber  and  may  be  spatially 
undersampled.  The  undersampling  is  performed  by  spacing  the 
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Figure 


3.9  Wavenumber  spectrum. 


phase  centers  of  the  subarrays  several  sensors  apart.  This  is 
equivalent  to  convolving  the  full  array  frequency  domain  data  vector 
with  the  subarray  conventional  beamforming  vector  (the  impulse 
response  of  the  FIR  filter)  and  then  undersampling  the  convolved 
output.  The  transfer  function  or  beampattern  of  a  subarray  is  shown 
on  the  wavenumber  plot  at  the  design  frequency  of  the  array  in 
Figure  3.10. 


The  true  spectrum  of  the  energy  impinging  on  the  array  is  thus 
passed  through  the  filter  defined  by  the  subarray  conventional 
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beamforming  vector.  The  filter  applies  an  amplitude  and  linear 
phase  change  to  the  spectrum.  Several  wavenumbers  centered  about 
the  presteer  direction  (the  MRA  of  the  subarray  beam)  are 
evaluated.  If  a  signal  has  a  wavenumber  that  is  not  on  the  MRA  of 
the  subarray  beam  it  will  be  attenuated  slightly.  The  wavenumber 


Subarray 

beampattem 


Frequency  (Hz) 

Figure  3.10  Beampattem  applied  to  wavenumber  spectrum. 


that  are  evaluated  are  in  the  mainlobe  of  the  subarray  beampattem 
and  receive  a  limited  amount  of  attenuation.  It  is  desired  to  avoid 
altering  the  signals  impinging  on  the  array  significantly.  This 
requires  that  several  adjacent  bands  of  wavenumber  be  filtered  (i.e. 
several  subarray  beams  or  presteer  directions)  and  processed  such 
that  the  total  wavenumber  spectrum  of  concern  is  evaluated  with  an 
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acceptable  amount  of  degradation.  The  attenuation  due  to  the 
subarray  beamforming  is  known  as  scalloping  loss. 

The  undersampling  of  the  spatial  spectrum  causes  the  aliasing 
previously  mentioned.  The  repeated  spectra  of  the  wavenumber 
spectrum  is  shifted  down  into  the  original  spectra.  This  occurs 
because  the  new  spatial  sampling  frequency  is  less  than  the  old  one 
due  to  the  greater  separation  of  the  data  samples  (phase  center  of 
the  subarrays).  Figure  3.11  shows  how  the  subarray  beampattern  of 
the  repeated  spectra  is  aliased  down  into  the  primary  spectra  on  the 
wavenumber  plot.  Figure  3.12  shows  two  adjacent  spectra 
beampatterns  and  their  exact  overlap.  Note  that  this  occurs  for  all 
adjacent  pairs  of  the  repeated  spectra. 


Aliasing 

spectrum 


Frequency  (Hz) 

Figure  3.11  Beampattern  aliasing  with  spectrum. 
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The  dimension  of  the  adaptive  beamforming  performed  on  the 
subarray  beam  outputs  is  desired  to  be  small  as  stated  in  the  first 
two  chapters  in  order  to  reduce  the  numerical  intensity  of  the 
algorithm  estimating  the  inverse  of  the  CSDM.  This  is  in  turn 
reflected  into  a  desire  for  as  few  subarrays  as  possible  that  will  span 
the  full  array  aperture.  This  calls  for  large  spacings  between 
subarrays  which  translates  into  a  very  small  spatial  sampling 
frequency.  This  in  turn  means  that  it  is  desired  to  have  the  repeated 
spectra  aliased  down  as  far  as  possible  into  each  adjacent  spectra. 


Aliasing 


-Ks  0  Ks 

Wavenumber 

Figure  3.12  Beampattern  aliasing. 


The  high  aliasing  is  not  desired  because  it  corrupts  the 
estimate  of  the  wavenumber  spectrum,  the  beam  output.  The 
compromise  between  the  two  conflicting  desires  was  suggested  by 
Dr.  Norman  Owsley  of  the  Naval  Underwater  Systems  Center  (NUSC) 
to  be  to  allow  the  aliased  spectrum  to  enter  an  adjactnt  one  such  that 
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the  two  mainlobes  overlap  at  the  level  of  the  first  sidelobe.  This 
condition  may  be  seen  in  Figure  3.12. 

The  question  arises  as  to  how  to  determine  what  subarray 
sizes,  spacings  and  shadings  meet  the  above  aliasing  criteria.  Due  to 
the  discrete  sensor  location  sampling  of  the  spatial  field  it  is  required 
that  there  be  an  integer  number  of  sensors  between  the  phase 
centers  of  adjacent  subarrays.  In  order  that  all  available  data  be 
used  it  is  also  required  that  the  subarrays  utilize  the  full  aperture  of 
the  array,  i.e.  that  the  first  and  last  sensors  are  included  respectively 
in  the  first  and  last  subarrays. 

Figure  3.13  depicts  the  integer  number  of  sensor  spacing  of  the 
subarray  phase  centers  and  the  utilization  of  the  full  aperture 
resulting  in  the  following  equation, 

M  =  ~ 1  ,  3-30 

JNa 

where  M  is  the  number  of  subarrays,  Ns  is  the  number  of  sensors  per 
subarray,  and  Na  is  the  number  of  sensors  between  the  phase  centers 
of  two  adjacent  subarrays.  Note  that  M,  N,  Ns,  and  Na  must  be 
integers.  The  number  of  sensors  in  the  full  array,  N,  is  given.  For  a 
particular  window  function  used  to  shade  the  subarrays  the  spacing 
of  the  subarrays,  Na,  for  a  particular  subarray  size,  Ns,  must  be  tested 
to  see  if  it  is  an  integer.  If  it  is  then  the  number  of  subarrays 
required  to  span  the  entire  aperture  of  the  array,  M,  is  tested  to  see 
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if  it  is  an  integer.  The  resulting  valid  combinations  will  satisfy  the 
constraints  and  conditions  previously  described. 


1  2  3 

Na  Sensors 
between  subarray 
phase  centers 

Figure  3.13  Subarray  spacing  example. 

Looking  at  the  exponent  of  the  conventional  beamforming 
vector  used  to  estimate  the  wavenumber  spectrum,  as  done  in 
Section  3.3. 1.1,  but  evaluating  it  at  the  spatial  frequency  ks  with  half 
wavelength  sensor  spacing  (d  =  Xs/2)  yields, 

bi(ks)  =  ks  y  i  i  =  0,...,N-1 

_  iSs.  Na^s  . 

"Na  2  1 


where  the  new  spatial  sampling  period  is  seen  to  be  Xu  =  NaXs/2  and 
the  corresponding  sampling  wavenumber  rate  is  ku  =  ks/Na.  As  seen 
in  Figure  3.12  the  new  sampling  rate,  ku,  is  exactly  two  times  the 
point  at  which  the  two  beampatterns  overlap  in  the  mainlobe. 


N  sensors  in 
the  full 
array 
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The  test  for  the  validity  of  a  particular  subarray  size  and 
shading  will  be  the  condition  that  the  ratio  of  the  original 
wavenumber  sample  rate  to  the  new  sample  rate  is  an  integer, 

Na  =  ri  =  INTEGER.  3-31 

Ku 

The  ratio  may  be  determined  by  computing  the  beampattern 
and  finding  the  wavenumber  in  the  mainlobe  that  is  at  the  same 
level  as  that  of  the  first  sidelobe.  From  this  the  new  wavenumber 
sampling  rate  may  be  determined  and  the  subarray  spacing 
determined. 

The  beampattern  must  also  be  evaluated  to  find  the 
wavenumber  that  is  at  the  limit  of  the  allowable  scalloping  loss. 
Once  this  is  determined  the  number  of  required  subarray  beams  and 
their  particular  presteering  directions  may  be  determined.  Figure 
3.14  shows  the  mainlobe  of  the  subarray  beampattern  and  the  region 
that  will  be  evaluated  where  the  scalloping  loss  is  acceptable,  [-ksi, 
ksll. 


The  corresponding  wavenumber  that  will  be  evaluated  on  each 
subarray  beam  will  be  less  than  the  band  shown  such  that  there  are 
an  integer  number  of  beams  that  span  the  full  wavenumber  space  as 
shown  in  Figure  3.15. 

The  constraint  is  seen  to  be  in  the  form  of, 


I  I  I  I  I  I 

-Ks/2  Ks/2 

Wavenumber 

Figure  3.15  Spacing  of  subarray  beams  in  wavenumber. 

The  band  of  wavenumber  (kc  wide)  that  is  evaluated  for  each 
subarray  beam  steering  direction  does  not  have  to  terminate  as 
shown  in  Figure  3.14  where  the  scalloping  loss  reaches  the  maximum 
allowable.  It  is  desired  to  have  the  same  width  of  wavenumber 
evaluated  for  all  the  subarray  beam  presteer  directions  and  to  have 
the  directions  be  equally  spaced  in  wavenumber. 
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The  algorithm  described  in  Appendix  C  was  used  to  determine 
the  valid  subarray  sizes  for  the  beampattern  and  beam  response 
evaluations  of  Chapter  4. 

Subarray  preprocessing  is  similar  to  a  vernier  frequency 
estimator  discussed  by  Nuttall  in  [15].  For  the  subarray  beamformer 
the  second  or  vernier  stage  of  processing  is  an  adaptive  spatial 
spectrum  estimator  instead  of  a  conventional  Fourier  method. 

3.3.2  Beam  Space  Preprocessing 

Beam  space  preprocessing  was  first  introduced  by  Gray  in  [3]. 
It  is  also  known  as  adaptive  beam  interpolation  (ABI).  The  sensor 
data  is  first  converted  into  N  (the  number  of  sensors)  conventional 
beam  outputs  performed  using  uniform  array  shading.  The  MRA’s  of 
the  N  conventional  beams  are  equally  spaced  in  wavenumber.  This 
insures  that  for  half  wavelength  sensor  spacing  the  interpolation 
beams  are  orthogonal  (diHdj  =  8jj). 

Several  adjacent  beams  are  then  used  to  interpolate  in  high 
resolution  the  area  at  the  center  of  the  fan  of  beams.  The 
interpolation  is  performed  with  a  distortionless  response  constraint 
as  was  developed  in  Section  3.3.1.  The  fan  of  beams  "wraps  around" 
when  the  angles  of  +90  or  -90  degrees  from  broadside  are  reached. 
This  is  allowed  because  of  the  repeated  spectra  generated  in  the 
spatial  sampling  process  as  discussed  in  Section  3. 3. 1.1. 
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The  mathematical  development  of  the  beam  space  algorithm  is 
covered  by  Gray  in  [3].  Reference  [4]  shows  how  the  beam  space 
algorithm  may  be  implemented  in  the  form  of  a  matrix  preprocessor. 
The  preprocessor  matrix  operating  on  the  array  sensor  data  is 
formed  from  the  steering  vectors  that  generate  the  conventional 
beams.  The  dimension  of  the  reduced  ABF  process  is  equal  to  the 
number  of  beams  used  in  the  interpolation.  The  number  of  presteer 
directions  is  equal  to  the  total  number  of  conventional  beams  utilized 
in  the  interpolation  process  which  is  equal  to  N,  the  number  of 
sensors. 

A  simpler  but  more  specific  implementation  of  the  beam  space 
preprocessor  could  take  the  form  of  first  the  computation  of  the  N 
conventional  beams  and  storage  of  the  complex  beam  outputs.  The 
beam  outputs  used  for  the  adaptive  algorithm  from  each  presteer 
direction  may  then  be  picked  at  will.  The  computational  savings  of 
this  implementation  over  the  matrix  preprocessing  is  realized  by  not 
having  to  recompute  conventional  beam  outputs  for  new  presteer 
directions. 

This  brings  to  bear  a  convenient  trait  of  the  beam  space 
preprocessor,  being  that  access  to  the  element  level  data  is  not 
required.  If  the  output  of  a  conventional  beamformer  is  the  only 
data  available  a  high  resolution  beam  response  estimation  may  still 
be  performed. 
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3.4  Enhanced  Minimum  Variance  Beamforming 

In  [1  &  10]  Owsley  introduces  the  enhanced  minimum  variance 

beamformer  (EMVDR).  At  low  interference  to  noise  ratios  the  MVDR 

adaptive  beamformer  performs  the  same  as  the  conventional 

beamformer.  This  can  be  seen  from  the  array  gain  comparison 

« 

shown  in  Section  4.1  and  [4].  Thus  the  ABF  achieves  most  of  its  gain 
against  the  stronger  interferers  impinging  on  the  array.  Owsley 
suggests  that  the  CSDM  be  approximated  by  the  components  due  to 
the  stronger,  or  dominant,  signals.  This  is  performed  by  separating 
the  frequency  data  vector  information  into  a  dominant  signal 
subspace  and  a  complementary  space  containing  all  other  signals. 
The  dominant  signal  subspace  is  then  scaled  or  enhanced. 

Mathematically  the  development  of  the  EMVDR  beamformer 
comes  from  the  eigenvalue-eigenvector  decomposition  of  the  CSDM. 
Recall  from  Section  3.1  the  CSDM  was  represented  as, 

R  =  D  S  DH  +  o2  IN  .  3  -5 

The  eigen-decomposition  may  be  represented  as, 

R  =  MAMH ,  3-33 

where  the  columns  of  M  are  the  orthonormal  eigenvectors  of  R  and  A 
is  a  diagonal  matrix  of  the  eigenvalues.  As  s^en  in  Brogan  [6]  and 
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derived  from  the  definition  of  R ,  the  following  are  some  of  the 
properties  of  the  eigenvalues  and  eigenvectors  of  R, 

-  The  eigenvectors  of  D  S  DH  are  also  eigenvectors  of  R 

-  If  \jfi  are  the  eigenvalues  of  D  S  DH  then 

Xi  =  Yi  +  a2  are  the  eigenvalues  of  R 

-  D  S  DH  is  positive  semidefinite  so  that  Yi  ^  0  for  all  i 

-  Since  R  =  RH,  then  MAMH  =  or  A  =  AH,  meaning  that 

the  eigenvalues  of  R  are  real 

-  If  o2  and  the  source  powers  in  the  diagonal  matrix  S  are 

greater  than  zero  then  the  eigenvalues  of  R  are  all 
positive  and  greater  than  zero  (R  is  positive  definite) 


Taking  the  D  dominant  eigenvalues  and  associated  eigenvectors 
and  separating  them  as  follows. 


E  =  N*PNH 


3-34 


where  N  is  an  N-by-D  matrix  whose  columns  are  the  D  eigenvectors 
of  D  S  DH  corresponding  to  the  D  largest  eigenvalues,  Yi*  and  ^  is  a 
diagonal  matrix  of  those  eigenvalues. 

The  exact  number  of  sources  does  not  need  to  be  known,  only 
the  dominant  or  very  strong  ones.  Many  high  resolution  source 
subspace  beamformers,  direction  finders  or  spectrum  estimators 
require  that  the  exact  number  of  sources  be  known  in  order  for  the 
algorithm  to  work.  There  has  been  much  work  in  the  area  of 
estimating  the  number  of  sources  impinging  on  an  array  or  the 
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number  of  sinusoids  present  in  a  time  domain  signal.  Most  of  the 
work  requires  that  the  background  noise  be  spatially  uncorrelated. 
Angularly  extended  noise  generates  eigenvalues  that  may  be 
confused  with  the  plane  wave  source  eigenvalues.  Thus  in  the 
presence  of  angularly  extended  noise  the  eigenvalue  methods  of 
estimating  the  number  of  sources  fall  apart.  The  EMVDR  algorithm  in 
contrast  only  requires  the  number  of  strong  plane  wave  sources.  It 
is  expected  that  this  problem  will  not  be  as  difficult  as  requiring 
knowledge  of  the  exact  number  of  sources. 

Once  the  dominant  source  subspace,  E,  is  determined  an 
estimated  and  enhanced  CSDM  is  formed  as, 

R(e)  =  eE  +  o2In 

=  eN'PNH  +  a2lN  3-35 

where  e  is  the  enhancement  factor. 

The  estimated  CSDM  is  used  in  an  MVDR  algorithm  operating  on 
the  full  array  frequency  data  vector.  The  beamforming  vector 
equation  is  found  by  setting  the  preprocessor  matrix  equal  to  an  N- 
by-N  identity  matrix  in  the  MVDR  development  of  Section  3.3.1 
yielding, 

NR-ld 

w  “  dHR*1d  ’ 
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Replacing  R1  with  the  approximated  CSDM,  R  (e)*1,  the  beamforming 
vector  becomes, 


NR(c)-*d 
dHR  (e)-id 


3-36 


Note  that  the  inverse  of  the  CSDM  is  required  and  may  be 
found  by  using  the  matrix  inversion  lemma  from  Brogan  [6], 


(UAVH  +  B)-i  =  B-i  -  B*lU(VHB*1U  +  A-i^VHB-i, 

as  follows, 

R  (e)-l  =  (eN'PNH  +  o2IN)-l 

=^[in-n(Id  +  7'p-1)-1nh] 

=  ^  [In  -  N  B(e)  N«] 

a2 

where  B(e)  =  (Id  + 

v 

=  diag(Pi) 

if  pi  =  - - TT. - 7  i=l,...,D. 

1  +o2/(e\j/i) 


Now  post  multiplying  the  inverse  CSDM  by  the  steering  vector  d. 


R  (e)-‘  d  =  ^[IN-NB(e)NH]d 
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=  "T[d  -  N  B(e)  NHd] 

1  £ 

=  3U  -  ZfrniniHd]  3-37 

0  i=l 

where  N  =  [ni  112  •••  no]  are  the  dominant  eigenvectors.  Placing  this 
back  into  the  equation  for  the  beamforming  vector, 

D 

N  [d  -  £pininiHd] 

•  | 

w= - 1-~ - .  3-38 

N  -  XPildHni|2 
i=l 

This  form  of  the  beamforming  vector  is  very  revealing  in  that 
it  shows  how  a  conventional  beam  is  formed  with  the  steering  vector, 
d.  Any  interferences,  defined  by  the  projection  of  d  onto  the 
dominant  eigenvectors  nj,  are  then  removed  by  steering  nulls  toward 
the  sources  contributing  to  the  information  in  the  dominant 
eigenvectors. 

A  considerable  advantage  of  the  EMVDR  technique  is  that  it 
requires  very  few  numerical  operations  for  implementation.  Kaveh 
and  Yang  introduced  a  very  fast  algorithm  for  the  estimation  of  the 
eigenvalues  and  eigenvectors  of  the  source  or  noise  subspaces  in  [7]. 
The  particular  algorithm  used  is  covered  in  Appendix  B. 


CHAPTER  4  -  COMPARISON 


4.0  Introduction 

This  section  will  compare  the  reduced  adaptive  dimension 
beamformers  with  conventional  and  full  array  adaptive 
beamformers.  The  metrics  used  for  comparison  are  the  array  gain 
improvement  (AGI)  for  a  single  source,  single  interference  and 
uncorrelated  background  noise  case,  the  beampatterns,  beam 
responses,  and  the  numerical  computation  intensities  associated  with 
the  realization  of  the  beamformers. 

4.1  Array  Gain  Improvement 

For  the  case  of  a  single  plane  wave  source  and  a  single  plane 
wave  or  point  interference  (PIN)  with  spatially  uncorrelated 
background  noise  the  array  gain  for  the  adaptive  beamformers  may 
be  found  analytically.  The  CSDM  for  the  above  case  according  to  the 
models  discussed  in  Section  3.1  is, 

R  =  SididiH  +  S2d2d2H  +  o2In  ,  4-1 

where  Si  and  S2  are  respectively  the  source  and  interference  powers, 
di  and  d2  are  their  associated  steering  vectors,  and  o2  is  the 
uncorrelated  noise  power. 
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The  array  gain  (AG),  as  described  in  Section  3.3,  is  the  ratio  of 
the  SNR  at  the  output  of  the  beamformer  to  the  SNR  at  a  single 
sensor.  The  single  sensor  or  beamformer  input  SNR  is  seen  to  be, 

SNR®,“S2+1<j2  '  4‘2 

The  beamformer  output  power  (beam  response)  from  Section 

3.2  is. 


r(0)  =  w(0)HRw(0)  . 

Inserting  the  above  definition  for  the  CSDM  for  the  single 
source,  single  interference  and  uncorrelated  background  noise  case 
and  dropping  the  dependence  on  the  angle  0, 

r  =  wH(SididiH  +  S2d2d2H  +  o2In)w 
=  SiwHdidiHw  +  S2WHd2d2Hw  +  ct2whInw 
=  SiwHdidiHw  +  S2WHd2d2Hw  +  a2wHw 
=  SilwHdil2  +  S2lwHd2l2  +  a2wHw 


The  beamformer  output  SNR  is  thus  seen  to  be, 

_  _  SilwHdil2 

SNR0UT-S2lwHd2|2  +  C2WHW  • 


4-3 


The  array  gain  is  then. 
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(S2  +  a2)lwHdil2 
^  ~  S2lwHd2l2  +  c2wHw 

(r  4-  N)lwH<iil2 
”rlwHd2l2  +  NwHw  ’ 


where  r  is  defined  as  the  ratio  of  the  interference  power  to  the  noise 
power  at  the  output  of  a  conventional  beamformer. 


r 


s2  _ns2 

<r2/N  ”  o2  • 
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The  array  gain  improvement,  described  by  Owsley  in  [4]  &  [10], 
is  defined  as  the  ratio  between  the  array  gain  for  a  particular 
beamformer  and  the  array  gain  for  a  conventional  beamformer  with 
uniform  shading.  The  beamforming  filter  vector,  w ,  for  the 
uniformly  shaded  conventional  beamformer  is  simply  the  steering 
vector  d .  This  beamformer  is  steered  to  look  at  the  source  of 
concern,  thus  w  =  d  i .  The  array  gain  is  then  seen  to  be, 

AGcBF=frTL’  4-6 


where  the  following  simplifications  were  applied, 

diHdi  =  N 
diHd2  =  Na 
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Id  iHd2l2  =  N^aa*  ■  N^L, 
a  is  in  general  complex  and  0  <,  L  £  1 . 

Note  that  if  there  is  no  interference  (S2  =  r  =  0)  the  array  gain 
for  the  conventional  beamformer  is  equal  to  N.  An  alternate  form  of 
the  array  gain  may  be  found  by  manipulating  the  equation  into  the 
form  of  a  constant,  N,  plus  a  second  term  dependent  on  the 
interference  to  noise  ratio  r, 

AGcbf  =  N  +  ,(‘+'rNLL).  4-7 


This  form  reveals  how  the  array  gain  is  affected  by  the  interference 
to  noise  ratio,  r. 


The  array  gain  for  each  of  the  adaptive  beamformers  may  be 
found  by  inserting  their  beamforming  vectors  into  the  above  generic 
AG  equation.  The  AGI  is  then  found  by  forming  the  ratio  between 
the  adaptive  beamformer  array  gain  and  the  array  gain  for  a 
conventional  beamformer. 


AGI  = 


AGabf 
AGcbf  * 


4-8 


An  exorbitant  amount  of  algebra  is  required  to  place  the 
equations  into  forms  that  are  meaningful.  Some  of  the  results  are 
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seen  in  [4]  and  encapsulated  in  Table  4.1.  The  EMVDR  beamformer 
AGI  was  derived  allowing  for  a  variable  enhancement  factor. 

The  parameters  associated  with  the  AGI  equations  are  as  follows. 

Beam  space: 

M  =  Number  of  interpolation  beams 
L  =  Average  sidelobe  level  of  interference  in  M  beams 

e  =  m"l 

Subarrays: 

L  =  Normalized  interference  level  in  full  array 
Ls  =  Normalized  interference  level  in  subarray 
L  s  =  Aliased  subarray  center  response 

EMVDR: 

P  =  Signal  to  interference  ratio  =  S1/S2 
2P 

Y  = - 

[1  -  PI  +  V  [1  -  P]2  +  4PL 

e  =  Dominant  source  subspace  enhancement  factor 
NS?  . 

r  =  -  interference  to  beam  output  noise  power 

U=  1  +  erP(l+y) 


It  is  desired  to  compare  the  four  adaptive  beamformers 
(MVDR,  subarrays,  beam  space,  and  EMVDR)  to  a  conventional 
beamformer  using  the  AGI  metric.  In  order  to  illustrate  the 
properties  of  the  equations  found  in  Table  4.1  three  separate  cases 
will  be  covered.  The  number  of  sensors  is  taken  to  be  N  =  32. 
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Table  4.1  Array  gain  improvement  (AGI)  for  full  array  MVDR  and 
beam  space,  subarray,  and  EMVDR  reduced  dimension  beamformers. 


j  +  r2L(l-L) 


gl&UL 
1  +  1+rE 


Beamformer  I  Array  Gain  Improvement 


Full  Array 
MVDR 


Beam  Space 
MVDR 


Subarray 
Beam  Space 
MVDR 


Enhanced***  _ [(a-L)+(a.lhLPn)F(l+»L) _ 
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interference  level  in  the  full  array,  L,  to  the  level  of  the  beampattem 
at  the  point  where  the  PIN  falls  in  the  mainlobe.  For  Figures  4.1  and 
4.2  the  PIN  falls  at  the  -3dB  point  in  the  beampattem  of  a 
conventional  beam  steered  at  the  source  (L  =  -3dB).  This  means  that 
if  a  conventional  beamformer  were  steered  to  the  direction  that  the 
source  was  arriving  from,  a  signal  coming  from  the  direction  of  the 
interference  would  only  be  attenuated  by  3dB.  For  strong 
interferences  this  is  not  very  much  attenuation. 


Figure  4.2  Array  gain  improvement  (AGI)  for  a  mainlobe  point 
interference  (PIN)  with  low  signal-to-interference  ratio. 


In  Figure  4.1  the  signal  to  interference  ratio  is  only  P  =  -3dB. 
This  means  that  the  signal  power  is  3dB  below  the  interference 
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power.  Figure  4.2  shows  the  case  where  the  signal  power  is  30  dB 
below  the  interference  power  (P  =  -30dB). 

In  Figure  4.3  the  PIN  is  placed  in  the  sidelobe  of  the 
conventional  beamformer  beampattem  such  that  it  receives  30dB  of 
attenuation  (L  =  -30dB).  The  signal  to  interference  ratio  is  P  =  -lOdB. 


CBF  MRA INTERFERENCE-TO-NOISE  RATIO  (dB) 

Figure  4.3  Array  gain  improvement  (AGI)  for  a  sidelobe  point 
interference  (PIN)  with  intermediate  signal-to-interference  ratio. 


The  AGI  for  the  EMVDR  beamformer  was  derived  for  the  single 
source,  single  interference  in  uncorrelated  background  noise  case. 
Only  one  dominant  eigenvector  and  eigenvalue  were  used  for  the 
estimate  of  the  CSDM.  Figures  4.1  and  4.3  show  the  EMVDR 
beamformer  reaching  a  plateau  at  high  interference  to  beam  output 
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noise  power,  r.  The  signal  to  interference  ratio  is  -3  and  -10  dB 
respectively  in  the  figures.  The  signal  is  thus  much  higher  than  the 
background  noise  level  for  high  r.  Because  of  the  strength  of  the 
signal  it  should  also  be  assigned  an  eigenvector.  This  would  improve 
the  performance  of  the  EMVDR  beamformer  to  make  it  equivalent  to 
that  of  the  full  array  MVDR  beamformer.  When  the  signal  power  is 
much  less  than  the  interference  power,  as  seen  in  Figure  4.2,  or  at 
low  interference  to  beam  output  noise  ratios  the  one  eigenvector  that 
is  used  to  estimate  the  CSDM  contains  enough  information  to  provide 
identical  performance  of  the  EMVDR  and  the  full  array  MVDR 
beamformers. 

Important  conclusions  to  draw  from  these  graphs  are  that  the 
reduced  adaptive  dimension  beamformers  perform  nearly  as  well, 
identically  or  better  than  the  full  array  adaptive  MVDR  beamformer. 
This  implies  that  the  performance  may  be  the  same  using  a  reduced 
adaptive  dimension  as  when  using  the  full  MVDR  beamformer.  This 
statement  must  be  clarified  to  include  the  condition  that  the  reduced 
adaptive  dimension  provide  enough  degrees  of  freedom  to  deal  with 
all  the  directional  interferences  present. 

Figures  4.1  and  4.2  show  how  the  adaptive  beamformers 
provide  substantial  improvement  over  the  CBF  when  the  interference 
to  noise  ratio  is  high  (r  >  lOdB).  When  the  interfering  signal  is  buried 
in  the  background  noise  (r  <  OdB)  the  adaptive  beamformers  provide 
no  improvement  over  a  CBF,  however  there  is  also  no  degradation. 
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This  implies  that  for  low  level  or  nonexistant  interferences  the 
adaptive  beamformers  perform  the  same  as  conventional  ones. 

4.2  Beampattems 

A  common  misconception  about  the  beampattems  of  high 
resolution  beamformers  concerns  their  mainlobe  width.  CBF 
beampattems  are  often  compared  by  their  mainlobe  width.  An 
adaptive  beamformer  that  provides  greater  resolution  does  not 
necessarily  have  a  narrower  mainlobe  than  its  CBF  counterpart.  The 
effect  of  the  adaptive  beamforming  is  in  the  placement  of  the  nulls 
that  occur  between  each  of  the  lobes  in  the  beampattern.  The  MVDR 
beamformer  places  nulls  such  that  the  energy  due  to  all  the 
interferences  is  optimally  cancelled.  In  order  to  illustrate  this  the 
example  of  Figure  4.1  is  explored.  The  interference  is  at  the  -3dB 
point  of  a  conventional  beam,  the  signal  to  interference  ratio  is  -3dB, 
and  the  interference  to  noise  ratio  is  15dB.  The  beamformers  are 
pointed  at  the  source  which  is  placed  at  the  broadside  angle  of  zero 
degrees.  The  interference  is  at  an  angle  of  1.58  degrees. 

The  even  figures  between  4.4  and  4.14  show  the  beampattems 
for  all  arrival  angles  for  six  different  beamformer  configurations. 
The  odd  figures  between  4.5  and  4.15  are  an  expanded  plot  of  the 
area  where  the  source  and  interference  are  located.  The  beamformer 
configurations  are  described  in  Table  4.2. 
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The  beampatterns  are  formed  assuming  that  perfect  knowledge 
of  the  CSDM  is  available.  This  is  the  case  of  infinite  time  averaging  or 
expected  performance  for  a  stationary  scenario. 


As  seen  in  Figures  4.4  &  4.5  the  conventional  beamformer  only 
attenuates  the  interfering  signal  by  3dB.  All  of  the  adaptive 
beamformers  steer  a  null  into  the  vicinity  of  the  interference  as 
shown  by  Figures  4.6  through  4.15.  The  gain  the  beampattern 
applies  to  the  source  angle  of  zero  degrees  is  always  unity.  This  is 
because  of  the  distortionless  response  constraints  imposed  by  all  of 
the  adaptive  algorithms. 


Table  4.2  Beamformer  configurations  for  Figures  4.4-4.15. 


EUUX£2 

Beamforminc 

PgasriBiiftn 

4.4  &  4.5 

CBF 

Uniform  shading 

4.6  &  4.7 

Full  array  MVDR 

4.8  *4.9 

Subarray 

#  Sensors  per  subarray  =  8 

#  Subarrays  =  7 

Subarray  shading  =  uniform 

4.10  &  4.11 

Beam  space 

#  Interpolation  beams  =  5 

4.12  &  4.13 

EMVDR 

Enhancement  =  1 

#  Eigenvectors  =  1 

4.14  &  4.13 

EMVDR 

Enhancement  =  1 

#  Eigenvectors  =  2 

The  beampattern  generated  by  the  subarray  preprocessing 
beamformer  shows  how  the  subarray  beampattern  is  imposed  on  the 
true  spectra  of  the  energy  impinging  on  the  array.  The  resulting 
spectrum  is  then  undersampled  causing  the  aliasing  as  seen  in  Figure 
4.8.  The  first  groups  of  three  lobes  to  either  side  of  the  main  lobes 
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(-15  to  -30  and  15  to  30  degrees)  are  the  most  adjacent  repeated 
spectra  that  are  aliased  down  into  the  original  spectrum.  Note  that 
the  beam  space  preprocessor  produced  a  beampattern,  Figures  4.10 
&  4.11,  almost  identical  to  that  of  the  full  array  MVDR  beamformer. 

The  EMVDR  beamformer  in  Figure  4.13  shows  the  null  slightly 
closer  to  the  source  than  the  other  adaptive  beamformers.  This  is 
expected  because  the  single  eigenvector  used  in  the  enhanced 
estimate  of  the  CSDM  that  forms  the  beamforming  vector  contains 
information  about  both  the  source  and  the  interference.  If  the  signal 
power  were  much  smaller  than  the  interference  (it  is  only  3dB 
smaller  in  this  example)  it  would  be  expected  that  the  information 
about  the  source  contained  in  the  dominant  eigenvector  would  be 
small  compared  to  the  information  about  the  interference.  Figures 
4.14  and  4.15  show  the  beampattern  generated  by  using  two 
eigenvectors.  Between  the  two  eigenvectors  all  of  the  information 
about  the  interference  is  found  and  it  is  seen  that  the  beampattern  is 
identical  to  that  of  the  full  array  MVDR  beamformer. 

4.3  Beam  Responses 

This  section  contains  a  comparison  of  the  responses  of  each 
beamformer  to  a  given  scenario.  Various  scenarios  are  presented  in 
order  to  illustrate  some  of  the  properties  of  the  CBF,  the  full  array 
MVDR  beamformer,  the  subarray  and  beam  space  matrix 
preprocessor  beamformers  and  the  enhanced  MVDR  beamformer. 


MAGNITUDE  (dB) 
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CONVENTIONAL  BEAMPORMER  BEAMPATTERN 


Figure  4.4  CBF  beampattem,  uniform  shading,  all  arrival  angles. 


CONVENTIONAL  BEAMPORMER  BEAMPATTERN 


ANGLE  OF  ARRIVAL  (DEG) 

Figure  4.5  CBF  beampattem,  uniform  shading,  selected  arrival 

angles. 
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Figure  4.7  MVDR  beampattern,  selected  arrival  angles. 


FULL  ARRAY  MVDRBEAMFDRMER  BEAMPATTERN 


FULL  ARRAY  MVDR  BEAMFORMER  BEAMPATTERN 
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ANGLE  OF  ARRIVAL  (DEG) 

Figure  4.6  MVDR  beampattern,  all  arrival  angles. 


MAGNITUDE  (dB) 


ANGLE  OF  ARRIVAL  (DEG) 

Figure  4.9  Subarray  beampattern,  selected  arrival  angles. 


(dfl) 


73 


BEAM  SPACE  BEAMPORMER  BEAMPATTERN 


-SO  -60  -40  -20  0  20  40  60  80 

ANGLE  OF  ARRIVAL  (DEO) 


Figure  4.10  Beam  space  beampattem,  all  arrival  angles. 


BEAM  SPACE  BEAMPORMER  BEAMPATTERN 


ANGLE  OF  ARRIVAL  (DEG) 

Figure  4.11  Beam  space  beampattem,  selected  arrival  angles. 
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Figure  4.12  EMVDR  beampattern,  one  eigenvector,  all  arrival 

angles. 


EMVDR  BEAMTORMER  BEAMPATTERN 


ANGLE  OF  ARRIVAL  (DEG) 


Figure  4.13  EMVDR  beampattern,  one  eigenvector,  selected  arrival 

angles. 


(HP) 
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EMVDR  BEAMPORMER  BEAMPATTERN 


Figure  4.14  EMVDR  beampattern,  two  eigenvectors,  all  arrival 

angles. 


EMVDR  0EAMFORMEK  BEAMPATTERN 


ANCLE  OF  ARRIVAL  (DEC) 

Figure  4.15  EMVDR  beampattern,  two  eigenvectors,  selected 

arrival  angles. 
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The  beamformers  are  driven  by  random  data  vectors  that  are 
generated  according  to  the  model  discussed  in  Section  3.1.  A  plane 
wave  source  impinging  on  the  line  array  has  a  complex  random 
amplitude  scaling  the  frequency  steering  vector  pointing  to  the 
direction  of  the  source.  Uncorrelated  complex  background  noise  is 
added  to  the  frequency  data  vector.  Thus  the  beamformer  input  is  a 
frequency  data  vector  consisting  of  the  sum  of  the  steering  vectors  of 
each  plane  wave  source  scaled  by  a  complex  random  number  and  a 
complex  random  noise  vector.  A  sixteen  sensor  array  with  one-half 
wavelength  spacing  is  beamformed  at  the  design  frequency  of  the 
array  to  look  in  5*N  =  80  beams  equally  spaced  in  wavenumber  or 
sin(arrival  angle)  space. 

The  graphical  output  is  in  the  form  of  an  arrival  angle  versus 
time  waterfall  of  beam  responses.  One  trace  in  the  time  waterfall 
corresponds  to  the  beam  response  of  the  beamformer  steered  to  each 
of  the  80  look  directions  for  a  particular  time.  The  time  increases 
upward  in  the  waterfall  where  the  first  (bottom)  trace  corresponds 
to  the  beam  response  computed  from  the  fifth  update  of  the 
exponentially  averaged  CSDM  or  the  eigen-structure.  The  last  (top) 
trace  is  the  beam  response  after  fifty  more  updates  to  the  CSDM  or 
eigen-structure.  This  display  structure  is  chosen  to  illustrate  how 
the  beamformers  perform  against  random  data.  Underneath  the 
waterfall  plot  each  beam  response  plot  is  overlaid  providing  a  form 
of  an  averaged  beam  response. 
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Each  beam  response  trace  is  compressed  to  provide  greater 
visual  display  capability.  The  compression  formula  is, 

log  io(  1  +  mx/a) 

y-k  logio(l+m)  ’  4-9 

where  y  is  the  value  plotted,  x  is  the  beam  output  power,  m  =  8,  a  =  1 
and  k  =  10.  The  compression  allows  smaller  signals  to  come  through 
in  a  linear  fashion  (seen  by  approximating  the  compression  function 
by  a  first  order  Taylor  series).  Larger  signals  are  compressed  by  the 
log  function.  The  end  result  is  that  low  level  signals  may  be  seen  in 
the  presence  of  higher  level  ones. 


The  deflection  or  detection  index,  as  described  by  Urick  in  [16], 
is  computed  for  each  plane  wave  source  in  the  beam  responses. 
Deflection  is  defined  as. 


where  ms  and  mn  are  respectively  the  mean  values  of  a  signal  and 
noise  beam  output  where  the  averaging  is  over  time  and  o2  is  the 
variance  of  the  noise  beam  output.  The  signal  beam  is  the  beam  with 
an  MRA  nearest  to  the  direction  of  arrival  of  the  source  of  concern 
and  the  noise  beam  is  placed  such  that  little  corruption  occurs  due  to 
plane  waves  entering  the  sidelobes  of  the  beampattern. 
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Table  4.3  Beam  response  case  scenario  descriptions. 


Figure 

Number 

Sourci 

Strength 

(dB) 

i  One 

Angle 

(degrees) 

Sourci 

Strength 

(4B) 

i  Two 

Angle 

(degrees) 

Noise 

Power 

(dB) 

Noise 

Bt*m 

(degrees) 

4.16 

0 

0 

hbhi 

-3.0 

-85.0 

4.17 

mmmmm 

0 

-10.0 

15.0 

-3.0 

-85.0 

4.18 

0 

0 

-10.0 

45.0 

-3.0 

-85.0 

4.19 

HHHilHH 

-3.0 

-16 

4.20 

0 

-45.0 

0 

0-50.0 

-3.0 

-85.0 

Five  different  scenarios  are  used  to  compare  the  beamformers. 
The  specifics  about  each  case  are  described  in  Table  4.3. 

Figures  4.16  through  4.20  show  the  graphical  output  of  the 
beam  responses.  The  CBF  beam  response  is  always  shown  as  the 
baseline  for  performance  as  seen  in  the  (a)  position  in  each  figure. 
The  shading  of  the  CBF  beamformer  is  Chebyshev  with  30  dB 
sidelobes.  The  full  array  MVDR  beamformer  is  in  the  (b)  position, 
the  beam  space  beamformer  with  five  interpolation  beams  is  in  the 
(c)  position  and  the  subarray  beamformer  with  five  subarrays  of 
eight  sensors  each  and  0.5  dB  acceptable  scalloping  loss  is  in  the  (d) 
position.  The  subarray  beams  have  Chebyshev  shading  with  30dB 
sidelobes.  The  (e)  and  (f)  positions  have  the  EMVDR  beamformer 
beam  responses.  The  (e)  position  is  used  when  one  dominant  mode 
or  eigenvector  is  used  in  the  beamformer  and  the  (f)  position  for  the 
case  of  two  dominant  modes.  The  enhancement  is  always  kept  at 
unity  to  keep  the  deflection  parameter  equivalent  to  the  other 
beamformers. 


Table  4.4  contains  the  deflection  parameters  for  the  beam 
responses  shown  in  Figures  4.!6  through  4.20.  The  scenarios  with 
two  plane  wave  sources  have  deflections  calculated  for  both  sources. 
The  source  number  shown  in  the  first  column  of  Table  4.4 
corresponds  to  the  source  number  shown  in  Table  4.3.  The  EMVDR 
beamformers  are  with  one  and  two  eigenvectors. 

Table  4.4  Plane  wave  source  deflections  (dB) 
for  beamformer  outputs. 

Figurc/Sourccl  CBF  MVDR  BEAM  SUB  EMVPR-1  EMVDR-2 

4.16 _ 29.20 _ 30.69 _ 23.79 _ 30,54 _ 32.17 _ - 

4.17  /  1 _ 27.82 _ 29.80 _ 27.21 _ 28.93 _ 30.39 _ 31.31 

4.17  /  2 _ 14.87 _ 11.50 _ 8i14 _ 11.00 _ 10.64 _ 11.67 

4.18  /  1 _ 27.61 _ 29.52 _ 26.41 _ 28.83 _ 30.26 _ 30.16 

4.18  !  2 _ 9^39 _ 11.81 _ TA6 _ 10.85 _ 11.07 _ 11.97 

4.19  /  1 _ 28.29 _ 30.05 _ 27.17 _ 28.90 _ 32.00 _ 32.20 

4.19  /  2 _ 25.39 _ 25.24 _ 22.14 _ 24.48 _ 24.93 _ 26.49 

4^20  18.97  27^93  25.48  27.00  28.64  25.42 

Fifty  time  samples  are  used  to  estimate  the  means  and 
variances  of  the  noise  and  signal  beams.  The  time  samples  used  are 
the  ones  shown  in  the  figures.  The  first  time  sample  shown  is 
generated  from  the  fifth  update  of  the  CSDM  or  eigenvector 
estimator.  This  is  to  remove  the  deflection  computation  from  most  of 
the  effects  of  the  initialization  of  the  CSDM  and  eigenvectors. 

Several  important  points  may  be  drawn  from  the  examples 
shown  in  Figures  4.16  through  4.20  and  are  encapsulated  in  the 
following  summary: 
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•  Almost  all  of  the  adaptive  beamformers  resolved  the  plane  wave 
sources  to  a  finer  bearing  than  the  conventional  beamformer.  The 
exception  is  when  an  eigenvector  was  not  attributed  to  a  plane  wave 
source  and  conventional  beamforming  was  performed  in  the  EMVDR 
beamformers.  In  this  case  the  resolution  is  equivalent  to  that  of  a 
uniformly  shaded  conventional  beamformer. 

•  A  low  level  source  is  more  clearly  seen  in  the  presence  of  a  high 
level  interference  for  all  of  the  adaptive  beamformers  (see  Figure 
4.17).  The  deflection  of  the  low  level  source  in  Figure  4.17  is 
greatest  for  the  case  with  two  eigenvectors.  This  is  expected  because 
of  the  high  resolution  obtained  by  the  EMVDR  beamformer  for 
signals  that  have  eigenvectors  assigned  to  them.  The  next  largest 
deflection  is  the  CBF.  This  is  not  expected  and  the  larger  deflection  is 
attributed  to  energy  entering  the  mainlobe  or  sidelobes  of  the 
conventional  beamformer,  steered  at  the  low  level  source,  from  the 
interference  at  zero  degrees.  This  will  cause  the  mean  of  the  signal 
beam  output  to  increase  and  correspondingly  an  increase  in  the 
deflection. 

•  A  very  important  aspect  of  the  EMVDR  beamformer  is  apparent  in 
Figures  4.19(e)  and  4.20(e)  where  the  high  level  signals  that  have 
eigenvector  assignment  are  very  finely  resolved.  All  other  angles 
have  conventional  beamforming  characteristics  with  the  added 
removal  of  the  interfering  signal. 
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•  Figure  4.20  shows  a  scenario  with  a  moving  source.  The  source 
changes  bearing  from  zero  degrees  to  fifty  degrees  in  fifty  time 
samples,  a  rate  of  one  degree  per  time  sample.  The  CBF,  MVDR, 
subarray  and  beam  space  beamformers  track  the  source  movement 
well  as  seen  in  Figure  4.20(a)-(d).  The  EMVDR  beamformer,  in  Figure 
4.20(e),  with  one  eigenvector  (assigned  to  the  stationary  source) 
allows  conventional  beamforming  of  the  moving  source  and  adequate 
tracking.  When  an  eigenvector  is  assigned  to  the  moving  source,  as 
seen  in  Figure  4.20(f),  the  eigenvector  estimation  algorithm  would 
not  track  the  changing  eigenvector.  This  may  be  due  to  an  incorrect 
averaging  constant  that  controls  the  amount  of  new  data  that  is  used 
by  the  eigenvector  estimator.  A  detailed  analysis  of  the  eigenvector 
estimation  algorithm  is  required  before  any  firm  conclusions  may  be 
drawn.  Note  that  this  questionable  result  does  not  detract  from  the 
EMVDR  beamforming  algorithm  for  other  methods  of  obtaining  the 
eigen-structure  of  the  CSDM  are  available. 

•  In  [4]  the  beam  space  and  subarray  matrix  preprocessors  showed 
a  scalloping  that  occurred  on  the  infinite  time  averaged,  or  perfectly 
known  CSDM,  beam  response.  The  scalloping  is  due  to  the  subarray 
conventional  beam  formation  and  the  attenuation  incurred  on  any 
signals  not  arriving  on  the  MRA  of  the  subarray  beam.  In  the  beam 
space  algorithm  the  full  array  conventional  beams  in  the 
preprocessing  are  what  causes  the  scalloping.  Note  that  the 
scalloping  is  not  apparent  in  the  beam  responses  to  random  data. 


IUE 
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(d)  Subarray  (e)  EMVDR,  D=1  (f)  EMVDR,  D=2 

Figure  4.17  Beam  response  -  Single  source, 
mainlobe  interference,  SIR  =  -lOdB 
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(a)  CBF  (b)  Full  array  MVDR  (c)  Beam  space 


(d)  Subarray  (e)  EMVDR,  D=1  (f)  EMVDR,  D=2 

Figure  4.18  Beam  response  -  Single  source, 
sidelobe  interference,  SIR  =  -lOdB 
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(a)  CBF  (b)  Full  array  MVDR  (c)  Beam  space 


(d)  Subarray  (e)  EMVDR,  D=1  (f)  EMVDR,  D=2 

Figure  4.19  Beam  response  -  Single  source, 
mainlobe  interference,  SIR  =  -3dB 
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(d)  Subarray  (e)  EMVDR,  D=1  (£)  EMVDR,  D=2 

Figure  4.20  Beam  response  -  One  stationary  source, 

one  moving  source 


87 


4.4  Numerical  Intensity 

The  full  array  MVDR  beamformer  requires  on  the  order  of  N2 
numerical  operations  per  update  per  beam  as  shown  by  Owsley  in  [2] 
and  [4].  The  inversion  of  the  estimated  CSDM  is  the  driving  factor 
behind  the  numerical  intensity.  The  reduced  adaptive  dimension 
beamformers  all  alleviate  this  problem  to  varying  extents.  As  seen 
in  [4]  and  Table  4.5  the  numerical  complexity  is  a  function  of  the 
following  parameters: 

MVDR,  CBF,  and  all  others: 

N  -  Number  of  sensors  in  the  array 

Lb  -  The  total  number  of  beams  to  form,  nominally  =  5N 

Beam  space: 

M  -  The  number  of  interpolation  beams 
Subarrays: 

M  -  The  number  of  subarrays 

Ns  -  The  number  of  sensors  per  subarray 

Ls  -  The  number  of  subarray  beams  to  form 

EMVDR: 

D  -  The  number  of  dominant  modes  or  eigenvectors  used 

Table  4,5  Beamformer  complexity  expressed  in  the  number  of 
arithmetic  operations  per  update. 


Beamformer 

Numerical  Complexity 

CBF 

NLh 

MVDR 

N2  +  N2Lh  +  NLb 

Beam  space 

N2  +  2M?N  +  M2Lb  +  MLb 

Subarray 

NSMLS  +  2M2LS  +  M2Lb  +  MLh 

EMVDR 

DN  +  DNLb  +  NLb 
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The  numerical  gain  achieved  by  the  use  of  a  beamformer  other 
than  the  full  array  MVDR  beamformer  may  be  calculated  by  taking 
the  ratio  of  the  number  of  operations  per  update  for  the  full  array 
MVDR  beamformer  to  the  required  number  of  operations  per  update 
for  the  topic  beamformer.  For  each  of  the  beamformers  of  the  beam 
response  comparison  of  Section  4.3  the  numerical  gain  is  computed 
using  the  equations  found  in  Table  4.5  and  displayed  in  Table  4.6. 


From  this  comparison  it  is  seen  that  the  EMVDR  beamformer 
with  one  eigenvector  achieves  the  most  numerical  gain.  The  beam 
space  and  then  the  subarray  beamformers  follow  slightly  behind  it. 
The  adaptive  dimension  of  the  subarray  beamformer,  the  number  of 
subarrays,  was  taken  to  be  equal  to  that  of  the  beam  space 
algorithm,  the  number  of  interpolation  beams.  Fewer  subarrays  or 
fewer  interpolation  beams  will  provide  greater  numerical  gain  and 
varying  performance.  The  price  of  optimal  and  sub-optimal 
beamforming  is  seen  by  the  large  numerical  savings  achieved  by 
using  conventional  beamforming. 


Table  4.6  Numerical  gain  for 
beamformers  of  Section  4.3. 


Beamformer 

Gain 

CBF 

17.2 

MVDR 

1.0 

Beam  space 

6.4 

Subarray 

5.7 

EMVDR,  D=1 

8.5 

EMVDR,  D=2 

5.7 

CHAPTER  5  -  CONCLUSION 


S.l  Summary 

This  thesis  provides  a  general  introduction  to  the  concept  of 
conventional  and  adaptive  beamforming  in  Chapter  2  with  a  proposal 
for  reducing  the  dimension  of  the  adaptive  process  in  adaptive 
beamforming. 

The  theoretical  development  of  the  matrix  preprocessed 
minimum  variance  distortionless  response  (MVDR)  beamformer  is 
presented  in  Chapter  3  as  well  as  the  beam  space  and  subarray 
matrix  preprocessing  schemes.  The  enhanced  minimum  variance 
beamformer  is  developed  in  Section  3.4. 

A  comparison  of  the  three  topic  reduced  dimension  adaptive 
beamformers  with  the  full  element  space  MVDR  algorithm  and  a 
conventional  beamformer  is  shown  in  Chapter  4.  The  comparison  is 
accomplished  using  such  metrics  as  the  array  gain  improvement 
(AGI),  the  beamformer  beampatterns,  beam  responses  and  numerical 
complexities. 

The  AGI  and  beampattem  comparisons  have  both  shown  for 
the  single  source,  single  interference  in  uncorrelated  background 
noise  scenario  that  the  reduced  adaptive  dimension  beamformers 
provide  nearly  the  same  average  performance  as  the  full  array 
MVDR  beamformer.  The  beam  response  comparisons  have  confirmed 
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this  in  terms  of  the  deflection  parameter  and  by  visual  inspection  of 
the  graphical  output  of  the  beamformers  shown  in  Section  4.3. 

The  enhanced  minimum  variance  beamformer  (EMVDR) 
showed  equivalent  performance  (when  the  correct  number  of 
eigenvectors  is  used)  in  the  AGI  and  beampattern  comparisons.  In 
the  beam  response  comparisons  the  EMVDR  beamformer  consistantly 
outperformed  the  MVDR  beamformer  in  terms  of  deflection  and 
visual  inspection  of  the  graphical  output.  There  are  still  two  areas  of 
concern  associated  with  the  EMVDR  beamformer:  the  required 
estimation  of  the  number  of  plane  wave  interferences  and  further 
exploration  of  the  eigen-structure  estimator  used  and  other  potential 
estimators. 

The  subarray  and  beam  space  reduced  adaptive  dimension 
beamformers  performed  nearly  as  well  as  but  not  better  than  the 
MVDR  beamformer  in  terms  of  deflection  and  visual  inspection  of  the 
graphical  output.  The  beam  space  beamformer  has  an  advantage 
over  the  other  adaptive  beamformers  in  that  access  to  the  sensor 
level  data  is  not  required.  The  subarray  beamformer  may  prove 
more  tolerant  to  array  sensor  position  perturbations  than  the  other 
adaptive  beamformers. 
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5.2  Future  Work 

Several  areas  of  interesting  work  deserving  greater  depth  have 
been  highlighted  throughout  this  thesis.  The  following  paragraphs 
are  devoted  to  listing  these  areas. 

The  EMVDR  beamformer  proved  to  be  the  most  interesting  and 
also  the  most  unfinished  of  the  beamformers  evaluated.  One 
primary  concern  of  the  beamformer  is  the  number  of  dominant 
modes  required  to  achieve  the  best  performance.  The  AGI 
comparison  in  Section  4.1  shows  that  only  the  strong  plane  wave 
sources  need  to  be  attributed  eigenvectors.  If  fewer  eigenvectors  are 
used  the  performance  degrades  some  but  as  seen  in  the  beam 
response  comparisons  of  Section  4.3,  particularly  Figure  4.19(e),  it 
may  be  possible  to  start  out  with  fewer  eigenvectors  than  required 
and  increase  the  number  after  visual  evaluation  of  the  beam 
response.  The  noted  figure  clearly  shows  a  strong  plane  wave  source 
being  beamformed  conventionally  in  the  beam  response  plot.  Once 
this  is  attributed  an  eigenvector,  as  in  Figure  4.19(f),  the  beam 
response  does  not  exhibit  any  obviously  strong  plane  wave  sources 
with  conventional  beamforming.  It  may  be  possible  to  automate  this 
process  removing  the  human  interface.  The  estimate  of  the  number 
of  dominant  plane  wave  sources  must  be  constantly  updated  for  the 
scenario  may  be  changing. 

The  second  concern  of  the  EMVDR  beamformer  is  the  method  of 
estimating  the  eigen-structure  of  the  CSDM.  Ideally  the  eigenvectors 
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and  eigenvalues  of  an  exponentially  averaged  CSDM  are  desired. 
This  will  allow  equivalent  comparison  of  the  EMVDR  beamformer 
with  other  beamformers  utilizing  the  exponential  CSDM  averager. 
The  one  used  in  the  simulations  of  this  thesis,  as  described  by  Yang 
and  Kaveh  in  [7]  and  in  Appendix  B,  is  numerically  very  quick  and 
may  provide  the  desired  quality  but  requires  further  exploration. 

As  mentioned  in  the  conclusions  section,  the  subarray 
beamformer  may  be  more  tolerant  to  array  sensor  position 
perturbations  than  the  other  beamformers.  When  the  exact  position 
of  the  sensors  in  the  array  is  not  known  and  the  position  may  be 
estimated  for  sections  of  the  array  the  subarray  technique  may  be 
ideally  suited  to  the  beamforming  task  for  the  perturbed  array.  The 
subarray  size  may  be  matched  to  the  section  whose  location  is  known 
and  then  beamformed  conventionally.  Those  beam  outputs  are  then 
combined  adaptively  with  the  extra  information  about  the  section 
postitions. 

A  further  application  of  the  subarray  preprocessing  may  be  the 
use  of  the  EMVDR  beamformer  at  the  subarray  beam  output  stage 
instead  of  the  MVDR  process.  This  would  further  reduce  the 
numerical  intensity  but  would  require  a  good  algorithm  for 
estimating  the  number  of  dominant  plane  waves  in  each  subairay 
beam  output. 
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The  CSDM  being  estimated  in  the  matrix  preprocessed  reduced 
dimension  adaptive  beamformers  is  of  a  smaller  dimension  than  that 
of  the  full  array  MVDR  beamformer.  This  may  provide  a  reduced 
noise  variance  on  a  beam  output  due  to  the  uncertainty  of  the  CSDM 
estimation.  This  has  not  been  shown  in  this  thesis  and  is  still  a 
conjecture,  however  it  is  expected  that  the  estimation  of  the  smaller 
dimension  CSDM  should  require  a  shorter  time  constant  in  the 
exponential  estimator  than  the  estimation  of  the  full  CSDM.  Faster 
adaptation  or  lower  variance  beam  outputs  may  be  achievable 
because  of  the  shorter  time  constant.  Further  work  in  this  area 
should  explore  the  estimation  of  the  reduced  dimension  covariance. 

As  a  group  beamformers  are  subjected  to  queries  about  their 
performance  against  several  debilitating  scenarios  such  as  angularly 
extended  noise,  correlated  plane  wave  sources  and  scenarios  that 
change  over  time.  The  reduced  dimension  beamformers  may  contain 
a  facet  that  aids  in  the  battle  against  these  types  of  problems.  Only 
further  study  will  provide  answers  to  these  queries. 


APPENDIX  A  •  Cross  Spectral  Density  Matrix  Update  and 

Inversion 


UPDATE  EQUATION 

The  covariance  or  cross  spectral  density  (CSDM)  matrix  is  not 
known  a  priori  and  must  be  estimated.  It  is  defined  as  the  expected 
value  of  the  outer  product  of  the  frequency  data  vector  with  its 
complex  conjugate  transpose  (tranjugate).  In  order  to  reduce  noise 
associated  with  the  estimation  of  the  CSDM  some  sort  of  averaging  is 
required.  Block  averaging  estimates  the  CSDM  as  the  sample  mean  of 
the  outer  products  of  L  frequency  data  vectors, 

-  i  i 

R=r  X  *(i»0xH(i»0  .  A 1 

L  i=l 


This  method  requires  L  frequency  data  vectors  to  generate  one  CSDM 
estimate. 

An  alternative  method  is  to  perform  exponential  averaging  of 
the  stochastic  rank  one  approximation  to  the  CSDM  yielding, 

R(k)  =  (l-a)x(k,f)xH(k,f)  +  aR(k-l)  ,  A  2 

where  a  is  the  exponential  averaging  constant  and  varies  between 
zero  and  one.  This  algorithm  requires  an  initialization  value  for  the 
CSDM.  A  block  average,  as  shown  above,  the  identity  matrix,  or  an 
intelligent  estimate  based  on  the  dominant  plane  wave  sources  and 
the  structure  of  the  CSDM  may  be  used. 

The  minimum  variance  distortionless  response  (MVDR) 
beamforming  algorithm  requires  the  inverse  of  the  CSDM  multiplied 
by  a  known  vector.  A  matrix  inversion  requires  on  the  order  of  N3 
numerical  operations.  For  large  values  of  N  this  is  not  acceptable. 
Owsley  covers  in  [2]  and  [5]  the  use  of  the  Cholesky  square  root 
factorization  of  the  CSDM  in  the  update  equation,  A2,  above.  In  [2] 
Owsley  shows  that  the  triangular  Cholesky  square  root  factorization 
yields  an  order  of  N2  numerical  operation  algorithm  for  determining 
the  MVDR  beamforming  vector.  The  Cholesky  square  root 
factorization  of  the  CSDM  is  defined  as, 
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R(k)  =  C(k)CH(k),  A3 

where  C(k)  is  the  Cholesky  square  root  factor  of  the  estimated  CSDM 
at  time  k.  Brogan  shows  in  [6]  that  the  Cholesky  square  root 
factorization  of  a  symmetric  positive  definite  matrix  exists  and  is  a 
lower  triangular  matrix.  Applying  the  Cholesky  square  root 
definition  of  the  CSDM  to  the  equation  A2  above  yields, 

R(k)  =  (l-a)x(k,f)xH(k,f)  +  aC(k-l)CH(k-l).  A4 

Now  by  rewriting  the  above  equation  as  the  product  of  an  N-by-N+1 
matrix  and  its  tranjugate, 

R(k)  =  [  '/ 1  -a  x(k,f)  V“C(k-l)  ][  V  1-a  x(k,f)  V^C(k-l)  ] 

=  U(k)UH(k).  A  5 


Note  that  the  N+l-by-N  matrix  UH(k)  is  in  the  form  of  a  row  vector 
appended  to  the  top  of  an  upper  triangular  matrix. 


UH(k)  = 


Vl  a  xH(k,f) 

L  VaCH(k-l)  J 


A6 


If  the  diagonal  of  the  upper  triangular  matrix  were  "zeroed  out"  then 
the  bottom  row  would  be  a  zero  row  vector  and  may  be  dropped 
yielding  an  N-by-N  upper  triangular  matrix. 

The  diagonal  of  the  matrix  CH(k-l)  in  A6  may  be  "zeroed  out" 
by  the  use  of  a  set  of  Givens  rotations  (see  Golub  &  Van  Loan  [13]). 
The  Givens  rotation  is  an  orthogonal  matrix  premultiplier  that  has 
the  ability  to  zero  out  an  array  element.  The  Givens 
premultiplication  matrix  is  in  the  form  of  an  N+l-by-N+1  identity 
matrix  with  a  two-by-two  matrix  slid  along  the  diagonal.  The  two- 
by-two  matrix  is  chosen  to  annihilate  the  diagonal  elements  of  the 
CH(k-l)  matrix.  The  orthogonality  of  the  Givens  rotation  matrix 
allows  the  preservation  of  the  U(k)UH(k)  product  in  equation  A5. 
The  N  Givens  rotations  produce  the  following. 


R(k)  =  U(k)PiP2-PNPNH-P2HPlHUH(k)  A  7 

=  U(k)PPHUH(k), 

where, 

P  =  P1P2-PN 

PiPiH  =  IN+1  . 

It  is  desired  to  form  the  update  equation  of  the  CSDM  into  a 
Cholesky  square  root  factorization  form.  Once  all  the  diagonal 
elements  of  CH(k-l)  are  annihilated  the  updated  Cholesky  square 
root  factor  is  simply  the  first  N  rows  of  the  resulting  Givens  rotated 
product, 

CH(k)  =  First  N  rows  of  PHUH(k). 

In  order  to  calculate  PHUIi(k),  the  N  Givens  rotations  must  be 
performed.  Because  of  the  form  of  the  Givens  premultiplier  matrix, 
PiH,  only  two  rows  of  UH(k)  are  altered  at  a  time.  Thus  the  updated 
Cholesky  square  root  is  formed  by  recomputing  the  ith  and  i+lst 
rows  of  UH(k)  for  i=l  to  N.  Because  of  the  sparseness  of  the  upper 
triangular  matrix  in  UH(k)  this  procedure  requires  very  few 
numerical  operations. 

BEAMFORMING  VECTOR  GENERATION 

As  seen  in  Section  3.2  the  MVDR  beamformer  requires  that  the 
inverse  of  the  CSDM  post  multiplied  by  a  known  vector  be  computed 
in  order  to  find  the  beamforming  vector.  This  is  seen  by  equation 
3-25  to  be. 


w(9) _ NV1AHdW  -  3-25 

w  dH(8)ARp'1AHd(0) ' 

Note  that  the  CSDM  inverse  post  multiplied  by  the  vector 
formed  by  the  projection  of  the  steering  vector  onto  the  presteering 
matrix  appears  in  both  the  numerator  and  the  denominator.  Thus  in 
order  to  form  the  beamforming  vector  the  CSDM  inverse  post 
multiplied  by  AHd(6)  need  only  be  formed  and  placed  into  equation 
3-25. 
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Forming  the  matrix-vector  product  as, 
v  =  R^u 

where  u  is  known  and  v  is  the  desired  output.  Then  multiplying 
both  sides  by  the  CSDM, 

Rv  =  u. 

Now  replacing  the  CSDM  by  its  Cholesky  square  root  factorization  as 
seen  in  equation  A3, 

CCHV  =  U. 

Let 

r  =  CHv  A  8 


where  substitution  yields, 


Cr  =  u.  A9 

This  equation  (A9)  may  be  solved  for  the  unknown  vector,  r,  by  a 
forward  substitution.  Once  r  is  known,  a  backward  substitution  may 
be  used  to  find  the  desired  vector,  v,  through  equation  A8.  The 
forward  and  backward  substitution  algorithms  are  shown  by  Golub 
and  Van  Loan  in  [13].  They  are  used  for  the  solution  of  the  Ax  =  b 
equation  for  x  for  lower  and  upper  triangular  matrices  respectively. 

The  beamforming  vector  may  now  be  computed  with  an  overall 
order  of  N2  numerical  operations  to  accomplish  the  CSDM  exponential 
update  and  the  matrix  inversion  required  by  the  MVDR  beamformer. 


APPENDIX  B  -  Eigen-structure  Estimation 
INTRODUCTION 

The  enhanced  minimum  variance  (EMVDR)  beamformer 
requires  the  estimation  of  the  D  dominant  eigenvectors  and 
eigenvalues  of  the  cross  spectral  density  matrix  (CSDM).  This  may  be 
accomplished  by  a  brute  force  general  eigenvector  decomposition 
routine  performed  on  the  exponentially  averaged  CSDM  or  by 
exploiting  the  structure  of  the  CSDM  and  estimating  only  the 
eigenvectors  that  are  desired.  One  of  the  algorithms  that  exploits  the 
strucutre  of  the  CSDM,  as  described  in  Section  3.1,  is  shown  by  Yang 
and  Kaveh  in  [7]. 

The  algorithm  described  in  [7]  also  utilizes  the  fact  that  the 
eigenvectors  of  a  matrix  and  of  a  polynomial  function  of  the  matrix 
are  the  same.  The  eigenvalues  of  the  polynomial  function  of  the 
matrix  are  the  original  eigenvalues  operated  on  by  the  polynomial 
function.  The  particular  polynomial  function, 

h(R)  =  R-i, 

causes  the  eigenvectors  corresponding  to  the  D  largest  eigenvalues  of 
R  to  be  the  D  smallest  eigenvalues  of  h(R).  Yang  and  Kaveh  have 
shown  that  this  particular  estimation  procedure  performs  better  than 
the  estimation  of  the  D  dominant  eigenvectors  of  the  non-inverted 
CSDM 

INITIALIZATION 

Yang  and  Kaveh  suggest  initializing  the  estimate  of  the 
eigenvectors  with  the  first  D  columns  of  an  averaged  CSDM  where 
there  must  be  at  least  D  averages  to  provide  enough  linearly 
independent  vectors.  The  first  D  columns  are  orthonormalized  using 
the  Gramm-Schmidt  technique.  It  was  found  that  initialization  of  the 
eigenvectors  with  an  N-by-D  identity  matrix  required  many 
iterations  of  the  algorithm  for  convergance  near  the  true 
eigenvectors.  Alternatively,  as  suggested  by  James  Nuttall  of  NUSC, 
the  eigenvectors  may  be  initialized  to  the  array  steering  vectors 
pointing  to  the  direction  of  the  D  dominant  sources.  This  method 
proved  to  work  very  well. 
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ALGORITHM 

The  algorithm  required  to  perform  the  update  is  now 
presented  in  vector  notation. 

Iteration  k: 


Input: 


Step  1: 


Step  2: 


a  -  Averaging  constant 

x  -  Frequency  data  vector. 

Wi(k-l)  -  Eigenvector  estimates  at  time  k-1,  i=l, 

Xl  =  wiH(k-l)x 

2(l-a)  x  Xi* 
wp  =  (2a-l)w  i  (k-1)  + - - - 


...,  D 


wi(k)  = 


/WpHWp 


Do  j  =  2  to  D, 

X  =  X  -  Xj.iWj.i(k) 

Xj  =  wjH(k-l)x 

2(1 -a)  x  Xj" 

wp  =  (2a-l)wj(k-l)  + - ^ - 


wj(k)  = 


End  do 


lwpHwp 


Output:  Wi(k)  -  Eigenvector  estimates  at  time  k,  i=l,  ...,  D 

Xi  -  Eigenvalue  estimates,  i=l,  ...,  D 


APPENDIX  C  -  Subarray  Sizing  Algorithm 

introduction 

The  size  and  spacing  of  a  subarray  is  constrained  to  allow  only 
a  specified  amount  of  aliasing  and  to  utilize  the  full  aperture  of  the 
array.  For  a  given  type  of  subarray  shading  each  possible  subarray 
size  must  be  evaluated  for  validity  in  the  above  two  concerns.  The 
aliasing  and  full  aperture  usage  constraints  are  elaborated  on  in 
Section  3.3.1. 

This  appendix  is  intended  to  introduce  the  algorithm  used  to 
compute  the  valid  array  size  and  spacing  combinations  when  the 
array  shading  and  the  number  of  full  array  sensors  is  given.  The 
scalloping  loss  limits,  as  described  in  Section  3. 3. 1.2,  are  found  as 
well  for  the  valid  subarray  size  and  spacing  combinations. 

DISCUSSION 

The  algorithm  requires  an  educated  searching  routine  that 
must  locate  the  maximum  point  of  the  first  sidelobe  of  the  subarray 
beampattem  generated  by  a  given  subarray  size  and  shading.  Once 
this  is  done  the  wavenumber  in  the  mainlobe  that  has  attenuation 
equal  to  that  of  the  first  sidelobe  maximum  is  determined.  From  this 
value  the  subarray  spacing  is  determined  and  tested  to  see  if  it  is  an 
integer.  If  it  is  an  integer  then  the  number  of  subarrays  required  to 
span  the  full  array  aperture  is  tested  for  an  integer  value 
characteristic.  If  both  the  integer  constraints  are  satisfied  the 
subarray  size  and  spacing  are  a  valid  combination.  The  beampattern 
is  then  searched  to  find  the  wavenumber  at  which  the  scalloping  loss 
limit  is  not  violated.  This  will  determine  the  number  of  subarray 
beams  required  to  estimate  the  full  arrival  angle  space  (-90  to  +90 
degrees). 

ALGORITHM 

The  algorithm  is  as  follows: 

Input:  N  -  Number  of  sensors  in  full  array 

wind(I)  -  Subarray  shading  function,  returns  shading 
for  I  elements  in  column  vector 
Seal  -  Scalloping  loss  limit  in  dB 

/*  Loop  through  ali  possible  subarray  sizes  (Ns  =  3  to  N-l)  */ 

Ncnt  =  0 
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Step  1:  If  (Ns<N-2)  Then 

Ns  =  Ns  +  1 

/*  Get  shadings  for  current  subarray  size  */ 
s  =  wind(Ns) 

/*  Compute  beampattem  via  FFT,  convert  to  dB  */ 
bp  =  101ogi0(abs(FFT(s,NFFT))2) 

/*  Step  through  beampattern  searching  for  scalloping 
loss  limit,  first  minimum  and  first  sidelobe 
maximum  */ 

/*  Initialize  for  search  */ 

Iscal  =  1 
Bmin  =  False 
Bold  =  bp(l) 

I  =  0 

Step  2:  1  =  1+1 

/*  Check  for  scalloping  loss  limit  */ 

If  ((bp(I)<Scal)  AND  (Iscal=l»  Then 
Iscal  =  Ind 
End  if 

/*  Check  for  first  minimum  */ 

If  ((bp(I)<BCUR)  AND  (Qmin  =  False))  Then 
Qmin  =  True 

Endif 

/*  Check  for  first  sidelobe  maximum  */ 

If  ((Bold>Bcur)  AND  (Qmin  =  True))  Then 
SLL  =  Bold 

Else 

Bold  -  Bqjr 
Goto  Step  2 

Endif 

/*  Find  the  aliasing  point  in  the  mainlobe  */ 

/*  Initialize  */ 

I  =  0 

Step  3:  1  =  1+1 

If  (bp(I)>SLL)  Then 
Goto  Step  3 

Endif 
Ialias  =  I 
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/*  Test  for  validity:  Integer  constraints  */ 
Step  4:  Na  =  Nfft  /  2Ialias 

M  =  1  +  (N  -  Ns)/Na 
If  (M  =  Integer)  Then 
Goto  Step  5 

Else 

Goto  Step  1 
End  if 


/*  Compute  output  parameters  */ 

Step  5:  Ncnt  -  Ncnt  +  1 

Ns(NCnt)  =  Ns 
M(Ncnt)  =  M 
Na(Ncnt)  =  Na 
L(Ncnt)  =  ceiI(NFFr/2IscAL) 
Goto  Step  1 


Else 

Stop 

Endif 

Output:  Ns 

M 
Na 
L 

Ncnt 


-  Vector  of  valid  subarray  sizes 

-  Vector  of  number  of  subarrays 

-  Vector  of  subarray  spacings 

-  Vector  of  number  of  subarray  beams 

-  Length  of  above  vectors,  number  of  valid 

combinations. 
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